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1.  Introduction 


This  section  will  provide  an  introduction  to  the  oxidation  of  textile  composites.  The  general 
challenge  of  characterizing  damage  due  to  environmental  conditions  will  be  discussed,  in 
particular  the  effect  of  a  high-temperature  oxidizing  environment.  The  development  of  a 
multiscale/multiphysics  finite  element  program  will  also  be  discussed  which  will  provide  the 
user  with  tailored  tools  for  textile  composite  analysis  that  do  not  exist  in  commercial  software.  In 
particular,  the  challenge  of  simulation  of  oxidation  in  textile  composites  and  the  coupling  of 
oxidation  and  mechanical  damage  will  be  addressed.  A  proposed  framework  for  the  coupling  is 
presented. 

1.1.  Damage  due  to  environmental  conditions 

Woven  composite  structures  are  expected  to  undergo  a  range  of  hygrothermal  and  oxidizing 
environmental  conditions  during  their  service  life.  Environmentally  induced  degradation  of 
textile  composites  has  been  examined  experimentally.  However,  the  characterization  is  typically 
macroscopic.  For  example,  Luan  et  al.  [1]  studied  the  corrosion  of  a  C-SiC  composite  with  SiC 
coating  (SiC-C/SiC)  under  a  low  frequency  cyclic  stress  in  various  gas  atmospheres  of  oxygen, 
water  vapor,  and  sodium  sulfate  vapor  at  temperatures  from  1000  to  1300  °C.  A  model  for  the 
cyclic  stress  corrosion  mechanism  of  the  composite  was  proposed  from  the  experimental  study 
and  an  equation  to  predict  the  lifetime  of  the  composite  under  cyclic  stress  conditions  was 
derived  from  the  model.  Hale  [2]  characterized  the  strength  reduction  of  three  GRP  composite 
materials  as  a  function  of  temperature  and  testing  environment  (sea  water  and  crude  oil 
condensate).  In  neither  case  were  the  microscopic  damage  mechanisms  considered. 

Haque  and  Rahman  [3]  investigated  the  damage  development  in  woven  ceramic  matrix 
composites  under  tensile  and  cyclic  loading  at  elevated  temperatures.  The  tensile  strength  of 
SiC/SiNC  woven  composites  was  found  to  increase  with  increased  temperatures  up  to  1000°C. 
Elevated  temperature  was  found  to  have  a  remarkable  effect  on  the  fatigue  strength.  At  700°C, 
the  fatigue  strength  was  approximately  50  percent  of  the  ultimate  strength,  while  at  1000°C  it 
was  found  to  be  less  than  20  percent  of  the  ultimate  strength.  They  developed  rate  equations  for 
modulus  degradation  and  life  prediction  under  fatigue  loading  at  room  and  elevated  temperatures 
which  fitted  well  with  the  experimental  results.  In  some  cases,  the  success  of  the  application 
itself  depends  on  the  ability  of  the  composite  to  withstand  environmental  conditions.  For 
example,  cryogenic  propellant  tanks  fabricated  using  composites  need  to  be  able  to  avoid  leakage 
of  the  propellant  through  the  micro  cracks  in  the  composite  material.  Peddiraju  et  al  [4] 
simulated  the  leakage  of  gaseous  hydrogen  through  the  thickness  of  a  damaged  composite 
laminate  and  predicted  the  leakage  rate  at  room  and  cryogenic  temperatures. 

Polymer  matrix  composites  absorb  moisture  during  service.  This  can  lead  to  plasticization  of  the 
polymer  matrix,  alter  the  stress  state  and  degrade  the  fiber/matrix  interface  [5-7].  Due  to  this,  a 
good  understanding  of  the  moisture  absorption  and  desorption  behavior  is  important  for 
predicting  long-term  material  and  structural  performance.  Some  good  work  has  been  put  in  to 
investigating  the  thermal  conductivity  and  moisture  diffusion  behavior  of  polymer  matrix  woven 
composites.  Dasgupta  and  Agarwal  [8]  studied  the  thermal  conductivity  of  plain  weave 
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composites  using  a  homogenization  technique  and  were  able  to  achieve  very  good  agreement 
with  the  experimental  results.  Roy  et  al.  [9]  examined  the  effect  of  preexisting  matrix-cracks  on 
the  moisture  diffusion  behavior  of  a  5 -harness  satin  weave  composite  using  a  continuum  damage 
mechanics  approach  based  on  the  theory  of  irreversible  thermodynamics.  Li  et  al.  [10] 
investigated  the  moisture  diffusion  behavior  in  hybrid  woven  composite  laminates  using  a  simple 
ID  diffusion  model  to  simulate  the  effect  of  stacking  sequence  of  woven  plies  on  the  diffusion 
behavior.  Tang  et  al  [11]  studied  the  effect  of  tow  architecture  on  the  diffusion  behavior  in 
woven  composites.  This  helps  in  identifying  the  dominant  architectural  factors  that  affect  the 
diffusion  behavior  of  a  polymer  matrix  woven  composite.  Their  analysis  consisted  of  two  steps  - 
calculating  the  effective  diffusivity  of  the  fiber  tows  with  matrix  and  then  using  these  properties 
to  model  the  tow  with  the  corresponding  tow  architecture  in  the  woven  composite.  The  effective 
diffusivity  of  the  tows  was  calculated  using  3D  finite  element  micromechanics  [12].  The  effect 
of  irregular  fiber  distribution  was  taken  into  account  using  a  finite  element  based  ‘bi-zone’  model 
[13],  Simulations  of  moisture  diffusion  tests  for  a  3-ply  woven  hybrid  composite  were  performed 
and  found  to  be  in  close  agreement  with  experimental  results. 


1.2.  Effect  of  oxidation 

Oxidation  at  high  temperature  has  been  a  concern  for  a  long  time.  Of  course,  the  definition  of 
high  temperature  depends  on  the  material  system.  Carbon  fiber-reinforced  silicon  carbide 
composites  (C-SiC)  exhibit  excellent  mechanical  properties  at  temperatures  below  1650°C  and 
have  been  designed  and  developed  for  high-temperature  applications  such  as  the  high  thrust-to- 
weight  ratio  turbine  engines  and  reentry  thermal  protection  for  spacecraft.  However,  the 
mismatch  in  thermal  expansion  coefficients  between  the  carbon  fiber  and  the  SiC  matrix  induce 
matrix  and  seal  coating  microcracking  during  cooling  from  the  processing  temperature  [14], 
These  cracks  allow  for  oxygen  to  leak  in  and  react  with  the  carbon  fibers  at  temperatures  above 
400  °C  [15-17].  This  oxidation  in  turn  will  degrade  the  mechanical  properties  of  the  composite. 
Luan  et  al  [18]  examined  C-SiC  composites  being  oxidized  or  corroded  in  various  gas 
atmospheres  and  found  that  oxygen  was  the  major  factor  degrading  the  composite  under 
conditions  with  cyclic  stresses.  They  proposed  a  model  for  the  cyclic  stress  corrosion  mechanism 
from  the  experimental  results  as  well  as  an  equation  to  predict  the  lifetime  of  the  composite. 
Halbig  et  al  [19]  studied  oxidation  tests  of  C/SiC  composites  at  elevated  temperatures  and 
developed  a  model  that  simulates  the  diffusion  of  oxygen  into  a  matrix  crack  bridged  by  carbon 
fibers. 

Carbon-carbon  composites  are  designed  for  extremely  high  temperatures,  but  they  must  be 
protected  from  oxidation.  Various  researchers  have  studied  the  behavior  of  carbon-carbon  and 
proposed  schemes  for  oxidation  protection.  Ceramic  coatings  alone  do  not  provide  a 
comprehensive  barrier  against  oxidation  because  of  the  mismatch  between  the  coefficients  of 
thermal  expansion.  Due  to  this  mismatch,  cracks  form  in  the  coatings.  As  an  additional  form  of 
protection,  particulates  are  added  to  the  matrix  [20],  These  ‘inhibitor’  particulate  materials  are 
usually  boron,  boron  carbide  or  silicon  carbide.  Ochoa  and  Elliott  [21]  studied  oxidation  under 
isothermal,  cyclic  thermal,  and  thermo-mechanical  fatigue  conditions  for  inhibited  carbon- 
carbon  composites.  Mass  loss  and  material  property  degradation  assessment  was  undertaken  with 
subsequent  exploratory  nondestructive  testing  utilizing  dynamic  mechanical  analysis  (DMA)  and 
piezoelectric  ultrasonic  composite  oscillator  technique  (PUCOT)  techniques.  Degradation  in 
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shear  and  axial  moduli  were  measured  as  oxidation  progressed.  Lou  et  al.  [22]  examined  the 
effect  of  additives  on  the  mechanical  properties  of  oxidation-resistant  carbon/carbon  composites 
(C/C).  The  additives  used  in  their  test  included  silicon  carbide,  silicon  nitride,  and  metal  borides. 
These  additives  resulted  in  large  increases  in  flexural  modulus  and  strength.  Recently  Mazany  et 
al  [23]  filed  a  patent  on  oxidation  inhibition  of  carbon-carbon  composites.  Their  invention 
involves  two  steps:  (a)  contacting  the  carbon-carbon  composite  with  an  oxidation  inhibiting 
composition  composed  of  phosphoric  acid  or  an  acid  phosphate  salt,  at  least  one  aluminum  salt, 
and  at  least  one  additional  metal  salt  and  (b)  heating  the  carbon-carbon  composite  at  a 
temperature  sufficient  to  form  a  deposit  from  the  oxidation  inhibiting  composition  within  at  least 
some  of  the  penetrated  pores  of  the  carbon-carbon  composite. 

Schoeppner,  Pochiraju  and  Tandon  [24]  developed  a  multidisciplinary  approach  aimed  at 
predicting  the  performance  of  high-temperature  polymer  matrix  composites  (HTPMCs). 
HTPMCs  are  used  in  a  variety  of  aerospace  applications.  Pochiraju  et  al  have  performed  an 
extensive  review  of  the  state  of  the  art  in  predicting  thermo-oxidative  degradation  and 
performance  of  HTPMCs[25].  Unfortunately,  there  is  still  much  more  research  required  and  all 
the  underlying  mechanisms  for  the  predicting  the  behavior  of  these  materials  are  yet  to  be 
determined.  Characterizing  the  behavior  of  these  materials  is  not  trivial  [27-29]  and  very  time- 
consuming  and  in  some  cases,  reliable  methods  to  determine  certain  properties  do  not  yet  exist. 
Pochiraju  et  al  also  reviewed  the  effect  of  oxidation  and  aging  on  the  fibers  as  well  as  composite 
behavior.  Tandon  et  al  [24]  characterized  the  behavior  of  neat  PMR-15  resin  and  developed  a 
model  to  predict  the  thermo-oxidation  of  the  material.  Thermo-oxidative  aging  was  simulated 
with  a  diffusion  reaction  model  in  which  temperature,  oxygen  concentration  and  weight  loss 
effects  were  considered.  The  model  which  was  implemented  using  FEM  considered  diffusion, 
reaction  and  oxidation  of  the  resin  system.  The  model  developed  by  Pochiraju  et  al  [24-26]  is 
used  as  the  basis  for  the  oxidation  model  developed  in  this  work  and  is  discussed  in  detail  in 
Section  3.  They  also  used  the  finite  element  method  (FEM)  to  model  the  oxidation  behavior  in  a 
Graphite/PMR-15  composite  [26],  where  they  assumed  the  fiber  did  not  oxidize.  The  oxidation 
model  developed  by  Pochiraju  et  al  tends  to  be  very  computationally  intensive  and  most  of  their 
analyses  were  performed  at  the  fiber/matrix  scale.  Pochiraju  et  al  [30]  also  used  the  oxidation 
model  to  predict  the  evolution  of  stresses  and  deformation  in  HTPMCs  by  accounting  for 
thermo-oxidation  induced  shrinkage.  The  oxidation  model  and  the  non-linear  elastic  deformation 
analyses  are  coupled  using  information  obtained  by  experimental  observation  of  shrinkage  in 
neat  PMR-15  resin  under  aging  in  oxygen  and  argon. 

Roy  et  al  [31]  developed  a  multi-scale  model  based  on  micromechanics  and  continuum  damage 
mechanics  to  simulate  the  accelerated  fiber-matrix  debond  growth  in  a  unidirectional  HTPMC 
undergoing  oxidation.  The  model  was  used  to  predict  the  mechanical  behavior  of  a  laminate  in  a 
three-point  bending  test  incorporating  the  damage  caused  due  to  oxidation.  Wang  and  Chen  [32] 
developed  a  computation  micromechanics  approach  based  on  irreversible  thermodynamics  to 
obtain  constitutive  properties  of  HTPMCs  while  tracking  thermo-oxidative  reactions, 
microstructural  damage  and  thermo-mechanical  loading.  A  two-scale  homogenization  theory  is 
also  used  to  determine  macroscopic  behavior  of  these  composites.  They  also  stressed  the  need 
for  many  not  yet  available  thermal,  chemical,  mechanical  and  interphase  properties  and 
microstructural  parameters  in  order  to  accurately  predict  the  behavior  of  HTPMCs. 
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1.3.  Development  of  multiscale/multiphysics  finite  element  framework 

There  are  many  commercial  and  public  domain  software  packages  for  finite  element  analysis. 
However,  they  are  typically  not  designed  for  the  particular  challenges  one  will  face  when 
performing  detailed  3D  analysis  of  textile  composite  structures.  Textile  composites  have  multiple 
microstructural  scales  -  the  fiber/matrix  scale,  the  lamina  scale,  and  the  laminate  scale.  This 
complex  microstructure  of  textile  composites  makes  it  necessary  to  use  multiscale  analyses  in 
order  to  obtain  detailed  information  about  their  behavior.  Moreover  the  proposed  work  also 
studies  the  behavior  of  textile  composites  under  oxidizing  environments.  This  requires  a 
multiphysics  analysis  that  couples  damage  progression  analyses  with  oxidation  simulations. 
These  sorts  of  novel  analysis  methods  are  not  convenient  to  implement  in  commercial  FEA 
packages  due  to  the  restrictive  nature  of  these  software. 

A  finite  element  analysis  framework  called  ‘BETA’  has  been  developed,  which  is  a  successor  to 
the  existing  in-house  finite  element  code,  "ALPHA".  Alpha  has  been  used  for  static  linear  and 
nonlinear  thermo-mechanical  analysis  and  transient  diffusion  analysis  of  textiles.  Existing  tools 
formed  the  foundation  of  the  BETA  finite  element  framework.  Although  the  existing  code  was 
designed  to  be  quite  modular  and  extensible,  experience  has  shown  that  the  needs  of  those 
performing  detailed  analysis  of  textiles  is  quite  severe.  The  new  framework  has  several 
enhancements  over  the  existing  in-house  code  in  order  to  meet  the  needs  of  the  proposed  work. 
The  goal  is  to  design  a  robust  framework  that  can  be  enhanced  and  extended  in  the  years  to  come 
by  future  users  and  lives  beyond  the  term  of  this  research  work.  Towards  this  end,  the  software 
has  been  designed  using  an  object  oriented  philosophy.  This  incorporates  features  such  as 
inheritance,  polymorphism,  data  abstraction  and  encapsulation.  When  designed  properly,  this 
kind  of  programming  philosophy  makes  it  easier  and  more  convenient  to  maintain,  manage, 
modify,  extend  and  enhance  a  large  software  package. 

The  new  framework  makes  use  of  the  latest  hardware  improvements  such  as  multi-processor 
machines  which  are  very  common  now.  The  framework  will  also  be  portable  so  that  it  can  be 
used  on  both  the  Windows  as  well  as  UNIX/LINUX  environments.  The  developed  framework 
can  be  used  to  analyze  different  configurations  including  textile  composites  subjected  to  a  high 
temperature  oxidizing  environment.  The  framework  includes  tools  for  geometric  description, 
including  spatial  variation  of  material  properties,  mesh  development,  finite  element  solver,  and 
postprocessing.  It  also  provides  better  control  of  output  for  debugging  algorithms  and 
postprocessing  of  results.  A  more  detailed  description  of  the  framework  is  given  in  Ref. [33], 


1.4.  Simulation  of  oxidation  in  textile  composites 

Composite  structures  are  increasingly  being  used  for  high  temperature  applications  in  the 
aerospace  industry.  The  extreme  operating  environments  that  these  materials  are  subjected  to  can 
lead  to  chemical  degradation  including  oxidation.  It  is  important  to  understand  the  behavior  of 
these  materials  under  these  conditions  so  that  they  can  be  designed  better  and  provide  increased 
performance.  A  focal  problem  that  is  investigated  in  this  work  is  the  effect  of  oxidation  on  the 
mechanical  behavior  of  textile  composites.  This  will  involve  a  coupled  damage  progression 
analysis  that  accounts  for  the  effect  of  oxidation  on  the  engineering  properties  of  the  composite. 

A  precursor  to  the  coupled  damage  progression  analysis  is  the  oxidation  analysis  of  the 
composite  which  is  quite  complex  because  in  reality  the  fiber  and  matrix  both  have  their  own 
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response  to  high  temperature  oxidation  and  aging.  In  addition,  when  the  two  are  combined  to 
form  the  composite,  the  anisotropic  oxidative  response  is  even  more  complex  to  simulate 
because  of  the  fiber-matrix  microstructure.  Micro-cracks  and  damage  formed  at  the  interface 
between  the  fiber  and  matrix  affect  the  oxidative  response  of  the  composite.  The  task  of 
simulating  oxidation  of  textile  composites  requires  a  combination  of  various  strategies.  The 
underlying  oxidation  model  is  adopted  from  the  work  by  Pochiraju,  Schoeppner  and  Tandon[24- 
26]  who  have  used  their  model  to  simulate  the  oxidation  of  neat  PMR-15  resin  with  reasonable 
accuracy  compared  to  experimental  observations.  The  oxidation  behavior  is  represented  using  a 
set  of  transient  nonlinear  governing  equations  based  on  the  conservation  of  mass  equation  for 
diffusion.  The  oxidation  model  will  be  implemented  using  the  finite  element  framework  that  is 
developed  as  part  of  this  work.  The  finite  element  formulation  imposes  limitations  on  the 
element  size  and  the  time  step  size  which  make  the  simulation  computationally  intensive.  New 
strategies  have  been  developed  in  order  to  expedite  the  oxidation  analysis.  Moreover,  it  is  not 
practical  to  discretely  model  the  fibers  in  the  textile  composite.  Strategies  for  determining 
effective  oxidative  properties  have  been  developed  and  validated.  The  overall  goal  of  this  work 
has  been  to  develop  an  efficient  analysis  strategy  that  can  simulate  the  oxidation  behavior  in 
textile  composites  in  a  reasonable  time  frame. 


1.5.  Prediction  of  Damage  in  textile  composites  under  oxidation 

The  overall  goal  of  this  work  is  to  use  a  finite  element  framework  to  analyze  damage  progression 
in  textile  composites  due  to  the  combined  effects  of  oxidation  under  high  temperature  and 
mechanical  loads.  Determining  the  effect  of  high  temperature  oxidation  and  aging  on  the 
mechanical  behavior  of  composites  is  a  very  complex  and  challenging  problem.  There  are  a 
number  of  studies  in  the  literature  investigating  the  different  time-dependent  physical,  chemical 
and  mechanical  damage  mechanisms  [25,34-36]  as  well  as  experimental  characterization 
studies[37-42].  But  there  is  still  much  more  work  that  needs  to  be  done  in  order  to  reliably 
predict  the  composite  behavior  using  mechanistic  approaches.  The  damage  progression  analysis 
involves  performing  an  oxidation  analysis  that  simulates  the  diffusion  of  oxygen  into  the 
composite  and  tracks  how  much  the  material  has  oxidized.  The  simulation  of  oxidation  in  the 
textile  composite  is  one  of  the  goals  of  this  work  and  is  discussed  in  the  previous  section. 

The  analysis  is  a  one-way  coupled  problem  where  the  oxidation  is  assumed  to  affect  the 
mechanical  behavior  of  the  material  and  not  vice  versa.  A  constitutive  theory  is  used  to 
determine  the  amount  of  damage  in  terms  of  strength  or  stiffness  degradation  based  on  the 
oxidation  state  of  the  material  in  the  composite.  Figure  1.1  shows  a  schematic  that  illustrates  the 
coupled  analysis.  Both  the  oxidation  analysis  as  well  as  the  damage  progression  analysis  needs  to 
account  for  the  multiple  microstructural  scales  in  the  composite.  The  damage  will  not  affect  the 
oxidation  properties  in  the  current  implementation.  The  progressive  damage  analysis  will  track 
the  damage  state  in  the  composite  and  calculate  the  stress  state  in  the  composite  with  respect  to 
time  as  the  oxidation  progresses. 
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Damage  creates  new 
pathways  for  diffusion 


Degrade  mechanical 
properties 


Stress  analysis 


Tow  Architecture 


|fe  ;;  Fiber/Matrix 

Micromechanics 


Oxidation/Diffusion 
analysis 


Figure  1.1:  Schematic  illustrating  coupled  oxidation/thermo-mechanical 


The  coupled  analysis  model  is  used  to  investigate  a  focal  problem.  The  focal  problem  chosen  for 
this  work  is  a  Graphite/PMR-15  plain  weave  composite  laminate  that  is  loaded  uniaxially  to  a 
particular  strain  level  and  then  the  top  and  bottom  surfaces  are  exposed  to  oxygen  for  200  hours. 
The  laminate  in  the  simulation  is  assumed  to  be  at  288  C.  A  parametric  study  was  performed  to 
study  the  effect  of  the  number  of  plies  in  the  laminate  on  its  mechanical  behavior.  This  analysis 
model  lays  the  groundwork  for  fully  coupled  simulations  of  the  behavior  of  textile  composites 
under  combined  mechanical  loading  and  oxidation. 

In  summary,  this  research  will  focused  on  the  following: 

1)  Develop  a  coupled  analysis  model  using  the  finite  element  framework  that  will  couple  the 
oxidation  analysis  and  the  damage  progression  analysis. 

2)  Develop  a  constitutive  model  to  simulate  the  effect  of  oxidation  on  the  mechanical  properties 
of  the  tow  and  matrix. 

3)  Use  the  coupled  analysis  model  to  analyze  a  focal  problem 

a)  Simulate  mechanical  behavior  of  a  Graphite/PMR-15  plain  weave  laminate  under 
oxidation. 
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b)  Perform  a  parametric  study  on  the  effect  of  the  number  of  plies  on  the  mechanical  behavior 
of  the  configuration. 


1.6.  Summary 

This  section  discussed  damage  due  to  environmental  conditions,  in  particular  oxidation.  For  the 
oxidation  analysis  described  in  this  work,  a  user-developed  finite  element  framework  provides 
the  flexibility  and  freedom  to  implement  the  required  model.  This  section  also  provided  a 
literature  review  that  detailed  briefly  the  challenges  and  accomplishments  involved  in  predicting 
the  effect  of  environmental  conditions  on  the  behavior  of  composites.  The  overall  goal  of  this 
work  is  to  develop  a  multiscale/multiphysics  analysis  framework  that  can  be  used  to  study  the 
mechanical  behavior  of  textile  composites  under  oxidation.  Section  2  will  describe  the 
architecture  of  woven  composites,  and  the  particular  configurations  of  interest  in  this  study. 
Section  3  will  describe  the  theory  and  formulation  of  an  oxidation  model  that  is  used  with  FEM. 
Section  4  will  provide  the  basic  algorithm  for  implementing  the  oxidation  model  into  a  finite 
element  program.  Section  5  will  discuss  strategies  that  were  employed  to  expedite  analysis,  and 
Section  6  will  present  validation  of  these  strategies.  Section  7  presents  homogenized  oxidation 
properties  obtained  through  analysis.  The  infeasibility  of  discretely  modeling  composites  on  the 
fiber-matrix  scale  necessitates  the  use  of  homogenized  properties.  Section  8  presents  results  for 
oxidation  analysis  of  a  plain  weave  textile  composite.  Section  9  presents  results  for  an  oxidation- 
damage  analysis  with  a  one-way  coupling.  Finally,  Section  10  will  present  conclusions  and 
possible  future  work. 
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2.  Configurations  of  Interest 


Textile  composites  are  commonly  utilized  in  high  performance  structural  applications  such  as 
those  in  the  aerospace  industry.  Such  materials  exhibits  higher  strength-to-weight  ratios  than 
traditional  metallic  materials.  Furthermore,  the  use  of  the  textile  architecture  can  ease  in  the 
manufacture  and  layup  of  composite  structures  compared  to  that  of  one  dimensional  laminates. 
However,  this  advantage  is  accompanied  by  an  increased  complexity  in  the  response  of  the 
material  due  to  the  more  complex  architecture  compared  to  tape  laminates.  A  variety  of  weave 
types  exist  with  varying  complexities.  This  section  will  discuss  weave  types  in  general  and  focus 
on  details  of  a  plain  weave  composite.  Furthermore,  modeling  of  a  textile  composite  on  a  fiber- 
matrix  scale  is  impractical  from  a  computational  point  of  view.  Thus,  homogenized  properties  of 
a  tow  should  be  developed  such  that  the  woven  architecture  may  be  modeled  discretely  on  a  tow- 
matrix  scale  while  still  characterizing  the  behavior  of  the  actual  composite  material. 

2.1.  Weave  types 

Textile  composites  provide  a  more  effective  manufacturing  process  compared  to  tape  laminates 
and  a  range  of  weave  architectures  exist.  Figure  2.1  shows  a  variety  of  weave  types  that  may  be 
considered  for  a  textile  composite.  The  plain  weave  type  is  relatively  simple,  while  the  various 
harness  satin  weaves  and  twill  weaves  have  a  more  complex  architecture.  These  weaves  are 
considered  two-dimensional  in  the  sense  that  there  is  no  through-thickness  stitching  or  weaving 
of  tows.  However,  it  should  be  made  clear  that  two-dimensional  weave  types  do  have  a  three- 
dimensional  spatial  variation  of  architecture.  Due  to  the  relative  simplicity  of  the  architecture, 
the  plain  weave  composite  was  selected  for  this  study.  This  weave  type  also  has  a  much  smaller 
repeating  unit  cell  that  serves  to  significantly  reduce  the  size  of  the  problem  being  analyzed. 


5 -harness  satin 


4-harness  satin 
(Crow) 


Basket 


8 -harness  satin 


Figure  2.1:  Examples  of  weave  types  utilized  in  textile  composite  materials 
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2.2.  Description  of  the  plain  weave  unit  cell 

The  relative  simplicity  of  the  plain  weave  unit  cell  makes  it  an  ideal  candidate  for  studying  the 
effect  of  a  thermo-oxidative  environment  on  textile  composites.  This  section  describes  the 
physical  dimensions,  waviness  ratio  (WR),  and  tow  volume  fractions  for  the  plain  weave 
configuration.  Waviness  ratio  is  defined  herein  as  the  thickness  (h)  of  a  woven  mat  divided  by 
the  wavelength  (A,)  of  the  undulation  of  a  tow.  A  composite  laminate  can  be  constructed  by  the 
duplication  and  translation  of  a  basic  unit  cell  (see  Figure  2.2).  The  matrix  is  shown  transparent 
in  the  figures  to  reveal  the  tow  architecture.  The  term  M  will  represent  the  basic  unit  cell  for  a 
simply  stacked  plain  weave.  The  term  M  will  represent  the  mirror  of  M  about  x?  =  h.  Therefore, 
a  symmetrically  stacked  unit  cell  may  be  denoted  by  the  stacking  sequence  [M  / M],  In  all  cases 
the  plain  weave  tow  undulation  was  sinusoidal  and  the  cross-section  was  lenticular  with 
sinusoidal  shaped  boundaries.  For  the  configuration  considered  in  this  study  a  waviness  ratio  of 
1/3  and  fiber  volume  fraction  of  55.6%  were  utilized. 


Figure  2.2  Description  of  plain  weave  unit  cell 


2.3.  Description  of  the  fiber-matrix  unit  cell 

The  woven  architecture  for  a  plain  weave  composite  is  well  defined  in  the  previous  section. 
However,  in  reality  the  tows  are  composed  of  discrete  fibers  and  matrix  which  are  impractical  to 
model  in  the  woven  configuration.  Thus,  effective  properties  for  the  tows  are  determined  by  the 
consideration  of  a  fiber-matrix  unit  cell  which  may  be  duplicated  and  translated  to  produce  the 
configuration  of  a  tow.  Figure  2.3  illustrates  the  multi-scale  analysis  involved  when  considering 
textile  composites.  The  particular  fiber  volume  fraction  may  be  accounted  for  by  increasing  the 
size  of  the  fiber  relative  to  the  surrounding  matrix.  For  lower  fiber  volume  fractions  a  square  unit 
cell  is  sufficient.  However,  for  higher  fiber  volume  fraction  with  closer  packing  of  fibers  a 
hexagonal  unit  cell  is  necessary.  For  this  study  a  square  array  of  fibers  was  used  to  obtain 
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homogenized/effective  properties  of  the  tows.  The  process  of  obtaining  and  validating  effective 
properties  is  discussed  in  Section  7  of  this  report. 
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Figure  2.3:  Multi-scale  analysis  of  textile  composites 


2.4.  Summary 

This  section  has  presented  a  variety  of  weave  types  commonly  utilized  in  textile  composites.  A 
plain  weave  type  was  chosen  for  the  relative  simplicity  of  the  architecture  and  the  ability  to 
significantly  reduce  the  problem  size  at  hand  for  the  oxidation  analysis  of  the  configuration.  The 
parameters  used  to  describe  the  plain  weave  unit  cell  were  explained.  The  necessity  of 
considering  a  fiber-matrix  unit  cell  was  explained,  and  the  usage  of  these  configurations  are 
explained  later  in  this  report. 
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3.  Theory  and  Formulation 


This  section  discusses  the  theory  and  formulation  for  diffusion  and  oxidation  using  finite 
element  analysis.  The  general  issue  of  boundary  conditions  for  periodic  microstructures 
is  discussed  as  well  as  an  algorithm  for  addressing  damage  initiation  and  growth.  Then 
the  theory  and  finite  element  formulation  for  diffusion  and  oxidation  are  discussed  in 
detail.  An  approach  for  the  coupling  of  oxidation  to  mechanical  analysis  is  presented, 
and  constitutive  relationships  are  discussed. 

3.1.  Boundary  conditions  for  periodic  microstructures 

Periodic  configurations  can  be  analyzed  by  using  a  representative  volume  element 
(RVE)  or  unit  cell.  They  can  also  be  used  to  obtain  effective  properties  for  the  periodic 
configuration  or  microstructure.  The  unit  cell  is  a  region  within  the  microstructure  which 
can  be  used  to  generate  the  entire  microstructure  by  just  duplication  and  translation  of 
the  unit  cell.  Once  the  unit  cell  is  chosen  for  the  periodic  microstructure,  the  certain 
characteristics  can  be  determined  based  on  the  fact  that  each  of  the  unit  cells  will  behave 
in  the  same  manner.  For  elasticity,  the  periodic  conditions  state  that  the  displacement  of 
one  unit  cell  differ  from  the  other  unit  cells  only  by  a  constant  offset,  which  depends  on 
the  volume  averaged  displacement  gradients  [43-4].  Furthermore  the  strains  and  stresses 
are  identical  in  all  of  the  unit  cells.  This  can  be  expressed  as 


da)=Ui(Xa)  +  (  — )  dn 

a)  i\  a)  fi 

(3.1) 

<0* 

II 

+ 
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(*a+da)  =  <TU  OO 

(3.3) 

where  da  is  a  vector  of  periodicity  [43-44].  The  vector  of  periodicity  is  a  vector  from  a 
point  in  one  unit  cell  to  an  equivalent  point  in  an  adjacent  unit  cell. 

Additional  computational  savings  can  be  obtained  by  exploiting  symmetries  within  the 
unit  cell  [43].  The  periodic  boundary  conditions  and  symmetry  conditions  are  imposed 
by  using  multi-point  constraints  in  the  finite  element  analysis. 


3.2.  Damage  initiation  and  progression 

The  damage  progression  analysis  performed  in  this  work  is  based  on  a  continuum 
damage  strategy.  This  strategy  degrades  the  strength  or  stiffness  of  a  material  point  in 
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the  finite  element  mesh  based  on  the  failure  criteria  for  the  material  point.  The  damage 
analysis  differs  with  respect  to  linear  elastic  analysis  lies  in  how  the  constitutive 
relations  evolve  as  the  load  on  the  configuration  changes.  This  section  will  describe  the 
algorithm  for  the  damage  progression  analysis  used  in  this  work  followed  by  the  failure 
criteria  and  the  property  degradation  scheme  used. 

All  the  analyses  performed  in  this  work  assume  that  the  configuration  is  loaded  with  an 
increasing  strain  load.  Figure  3.1  shows  the  flowchart  for  the  algorithm  used  in  this 
work.  The  configuration  is  assumed  to  behave  linearly  until  damage  is  initiated.  The 
failure  criterion  can  be  used  to  determine  the  load  at  which  damage  initiates.  This  is 
done  by  applying  an  arbitrary  load  on  the  model  and  calculating  the  expressions  in  the 
failure  criterion.  Since  the  model  is  initially  linear  elastic  until  the  first  instance  of 
damage,  it  is  possible  to  calculate  the  initial  failure  load  by  scaling  the  load  to  the  point 
where  failure  is  just  triggered.  The  configuration  is  then  loaded  with  a  load  that  is  a 
small  fraction  larger  than  the  load  at  which  damage  initiated.  This  is  done  to  ensure  that 
round-off  errors  during  the  numerical  calculations  do  not  affect  the  simulation.  This 
ensures  that  the  load  is  increased  to  a  value  that  definitely  causes  new  damage  to  occur. 
This  point  in  the  algorithm  can  be  considered  the  beginning  of  a  new  load  step.  The 
displacement  field  corresponding  to  this  load  is  solved  for,  by  assuming  that  no  damage 
has  initiated  yet. 

The  displacement  field  is  used  to  apply  the  failure  criterion  at  all  the  Gauss  quadrature 
(or  integration)  points.  For  all  the  locations  that  damage  is  found,  the  mechanical 
properties  at  that  integration  point  are  degraded  based  on  the  property  degradation 
scheme.  The  model  is  solved  for  the  new  displacement  field  based  on  the  new  material 
properties  at  each  integration  point.  The  model  is  checked  again  for  damage  and  this 
procedure  is  repeated  until  no  new  damage  is  detected.  The  next  step  before  moving  on 
to  the  next  load  step  is  determining  the  load  for  the  next  load  step.  Since  we  have 
converged  to  a  damage  state  for  this  current  load  step,  the  configuration  can  be  likened 
to  a  new  linear  elastic  material  till  the  load  is  increased  and  new  damage  is  found.  Thus, 
just  as  the  load  for  initial  failure  was  determined,  the  load  value  for  the  next  occurrence 
of  new  damage  is  determined  using  the  failure  criteria.  In  this  manner,  the  load  is 
increased  and  the  simulation  proceeds  through  the  load  steps  until  a  specified  maximum 
strain  load  is  reached.  Throughout  this  process,  the  damage  state  is  recorded  and  new 
damage  is  tracked  as  the  load  on  the  configuration  is  increased.  Other  post-process  data 
such  as  the  volume  averaged  stresses  and  strains  are  also  recorded.  Figure  3.2  gives  a 
schematic  of  what  the  stress-strain  response  would  look  like  as  the  simulation 
progresses.  The  following  sections  describe  the  failure  criteria  and  the  property 
degradation  scheme  that  were  used  in  this  work. 
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Figure  3.2:  Schematic  of  stress-strain  response 


3.2.1.  Failure  criteria 

There  are  various  failure  criteria  such  as  the  maximum  strain  criterion  and  the  maximum 
stress  criterion.  Any  of  these  criteria  or  a  combination  of  these  criteria  can  be  used  in  the 
damage  progression  analysis.  For  the  analyses  in  this  research  work,  the  maximum  stress 
failure  criteria  are  used.  This  means  that  failure  has  occurred  when  any  of  the  stress 
components  in  the  material  coordinates  has  exceeded  its  corresponding  strength,  i.e. 
when  CTy  /  Sy  >  1  where  cr.  is  the  stress  component  in  the  material  coordinates  and  Stj 

is  the  corresponding  strength  for  cr. .  Section  9  gives  a  detailed  description  of  the  failure 
criteria  that  are  utilized  to  simulate  the  microscopic  damage  progression  in  this  work. 


3.2.2.  Property  degradation  scheme 

Typical  property  degradations  models  degrade  the  engineering  properties  whenever 
failure  is  detected.  Some  degradation  models  look  at  the  properties  (such  as  stress, 
strain)  at  the  center  of  the  element.  In  this  work,  the  failure  criteria  and  property 
degradation  scheme  is  applied  on  each  integration  point  of  all  the  elements  in  the  model. 
The  stresses  and  strains  at  any  material  point  in  the  material  coordinate  system  are 
related  by  Hooke’s  law  given  by  Eq.(3.4).  The  compliance  matrix  for  an  orthotropic 
material  is  given  by  Eq.(3.5). 
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(3.5) 


Let  Etj ,  Gij  and  vtj  be  the  original  extensional  moduli,  shear  moduli  and  Poisson’s 

ratio  respectively  and  Eu  ,  Gij  and  Vy  be  the  degraded  extensional  moduli,  shear 

moduli  and  Poisson’s  ratio  respectively.  Say,  cii ,  (i  =  1..9)are  the  degradation  parameters, 

which  specify  the  amount  of  degradation.  Then  a  typical  property  degradation  scheme 
will  look  like: 


£|  j  /l|  |  /  (- / 1  ,  £"22  ^  1  ^  ,  £33  £33  /  Cly 

G12  —  G12  /  fl4 ,  G23  —  G23  /  a5 ,  G13  —  G33  / 


12  ^12 


/ i  I  (I —  ,  t^23  t^23  /  1  ,  / 3|  i  /  1 1 1  /  1  / 


13 


33 


(3.6) 


For  example,  if  <2,  =8,  that  implies  that  the  £11  modulus  is  decreased  by  a  factor  of  8 

from  its  current  value  if  the  material  point  fails.  Note  that  in  this  general  framework,  the 
diagonal  as  well  as  non-diagonal  entries  of  the  compliance  matrix  can  be  affected 
independently.  The  specific  details  of  property  degradation  scheme  used  in  this  work 
including  the  degradation  factors  used  for  the  different  materials  will  be  given  in  Section 
9. 
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3.3.  Diffusion 


This  section  describes  the  governing  equations  and  the  finite  element  formulation  for 
simulating  the  diffusion  behavior  in  materials.  The  diffusion  behavior  through 
heterogeneous  materials  was  analyzed  in  this  work.  This  section  starts  by  listing  the 
equations  for  the  common  analysis  procedure.  This  is  followed  by  the  derivation  of  the 
weak  form  and  its  discretization  to  obtain  the  finite  element  formulation. 

3.3.1.  Conservation  laws 


The  conservation  of  mass  law  for  diffusion  yields  the  following  equation 


sc_  a j± 

dt  dxi 


0 


(3.7) 


where  C  is  the  concentration  of  diffusing  material  and  Jt  is  the  diffusion  flux. 


The  differential  equation  described  in  Eq.(3.7)  holds  for  a  material  point.  When  the 
material  being  analyzed  is  homogenous,  the  concentration  field  is  continuous  throughout 
the  domain  and  can  be  solved  without  any  modifications.  When  the  governing  equation 
is  applied  to  a  configuration  that  has  heterogeneous  regions  with  dissimilar  solids,  the 
concentration  is  generally  not  continuous  across  the  interface  between  the  different 
solids.  This  issue  of  discontinuous  concentrations  is  addressed  in  Ref. [12],  where  a 
thermodynamic  potential  is  introduced.  The  thermodynamic  potential  is  considered  to  be 
what  drives  the  flow  of  a  diffusing  material  through  another  material.  This  potential  is 
continuous  across  the  material  interface  and  the  concentration  is  then  defined  as  a 
function  of  the  thermodynamic  potential.  When  this  function  is  assumed  to  be  linear 
with  C=0  when  the  potential=0,  the  function  is  of  the  form 

C  =  aC  (3.8) 


where  C  is  the  thermodynamic  potential  and  a  is  a  material  property.  C  is  assumed  to 
have  a  range  from  0  to  1,  which  means  that  the  concentration  is  maximum  when  the 
potential  has  a  value  of  1.  That  determines  a  to  be  the  saturation  mass  concentration  of 
the  diffusing  material  in  the  solid,  denoted  by  C°° .  Therefore,  the  thermodynamic 
potential  is  the  concentration  in  the  material  normalized  by  the  saturation  concentration 
of  the  solid,  hereafter  referred  to  as  the  normalized  concentration, 

C  =  —  (3.9) 

Cx 

The  governing  equation  can  now  be  rewritten  as 
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(3.10) 


C"^£+A,I=0 

dt  dx. 


3.3.2.  Concentration  gradients 

Unlike  the  strain-displacement  relationship  in 


uses  simple  gradients  of  the  concentration 


solid  mechanics,  the  diffusion  analysis 


3.3.3.  Constitutive  relations 


The  relationship  between  flux  and  the  concentration  gradient  is  given  by  Fick’s  first  law, 


J:  =  -D: 


dC_ 

dx 


(3.11) 


where  Z>  is  the  2nd  order  diffusivity  tensor.  The  Latin  subscripts  i  and  j  denote  the 
coordinate  direction  and  range  from  1  to  3  for  a  three  dimensional  formulation. 

When  Eq.(3.11)  is  re-written  in  terms  of  the  normalized  concentration, 

dC 

J:  =  -C"Z).  — -  (3.12) 

‘  9  dx , 


3.3.4.  Boundary  conditions 

The  flux  boundary  conditions  are  given  by 

q  =  —niJi  onS  (3.13) 

And  the  normalized  concentration  boundary  conditions  are  given  by 

C  =  C  on  5  (3.14) 

where  C  is  the  specified  normalized  concentration  on  the  boundary  S 
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3.3.5. 


Weak  form 


The  equation  of  conservation  mass  is  multiplied  by  a  variation  of  the  normalized  concentration 
and  integrated  over  the  volume  to  obtain  the  weighted  integral  form, 


r  ^  — 

oo  dC  d 

\8C 

C  -  +  -  J: 

J 

y 

1 

Cl) 

1 _ 

IdV  =  0 


(3.15) 


where  SC  is  an  arbitrary  variation  of  the  normalized  concentration. 
Integration  by  parts  gives  the  weak  from, 


dSC 

dx. 


Ji  dV  +  \SCniJidS  =  0 


(3.16) 


Using  Eq.(3.12)  and  Eq.(3.13)  in  Eq.(3.16)  gives  the  basis  for  the  finite  element 
formulation, 


dC  dSC 
- 1 - 

dt  dxt 


CXD  — 

dx. 


W 


s 


(3.17) 


3.3.6.  Discretization  of  weak  form  and  time  integration 

Over  a  typical  finite  element,  the  normalized  concentration  is  approximated  by 

C(x,t)  =  Na(x)Ca(t)  (3.18) 


where  Na  are  the  interpolation  functions  and  Ca  are  the  nodal  normalized 

concentrations.  The  subscripts  with  Greek  letters  range  from  1  to  the  number  of 
interpolation  functions. 

After  discretizing  the  weak  form  using  Eq.(3.18)  and  SC  =  NaSCa ,  the  following  set  of 
equations  are  obtained, 


dCB 

A j  cm  C 

a  p  dt 


dN  dN„  - 
dx..  ,]  dx;  p 


W 


pjds 


(3.19) 


In  matrix  form  this  can  be  written  as 

M a[/lp  T  Kapqp  =  Fa  (3.20) 
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where 


V 

(3.21) 

rr  dN  dN„ 

Kafj  =  a  C  D  p  dV 

'  dxt  _ 

(3.22) 

K=jNJds 

s 

(3.23) 

R 

II 

S" 

a 

s 

II 

&  srii 

(3.24) 

Note  that  Eq.(3.19)  contains  a  time  derivative  term.  In  order  to  numerically  solve  this  set 
of  equations,  an  approximation  is  used  for  the  time  derivative  term  whereby  the  solution 
at  a  particular  instant  in  time  is  determined  based  on  the  solution  history.  The  following 
describes  this  ‘time-marching’  procedure  used  to  numerically  solve  Eq.(3.19).  For  the 
sake  of  convenience,  the  following  generally  uses  matrix  notation  and  dispenses  with  the 
indices.  Let  the  subscript  5  denote  the  solution  at  time  5  and  the  subscript  s+1  denote  the 

solution  at  time 5  +  At.  Using  Eq.(3.20),  the  following  equations 
consecutive  time  steps,  t=ts  and  t=ts+i. 

can  be  written  for  two 

Mqs+Ksqs-Fs=  0 

(3.25) 

MLi  Ks+lqs+i  ~  Fs+ 1  —  0 

(3.26) 

Using  the  alpha  family  of  approximations  [45]  gives 

(1  -  a)qs  +  aqs+l  =  Aq  /  At 

(3.27) 

Multiplying  Eq.(3.27)  by  A tM  gives 

(1  -  a)AtMqs  +  aAtMqs+]  =  MAq 

(3.28) 

Rearranging  the  terms  in  Eq.(3.28)  gives  an  expression  for  aAtMqs+1 

aAtMqs+l  =  M  Aq-(\ -  a)AtMqs  (3.29) 

Multiplying  Eq.(3.26)  throughout  by  aAt  gives 
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aAtMqs+1  +  aAt[Ks+1qs+1  - Fs+t  ]  =  0  (3.30) 

Substituting  Eq.(3.29)  into  Eq.(3.30)  gives, 

MAq  -  (1  -  a)AtMqs  +  aAt  [Ks+lqs+1  -Fv+1]  =  0  (3.31) 

And  substituting  for  Mqs  from  Eq.(3.25)  in  Eq.(3.31)  gives  the  governing  equations 

MAq  -  (1  -  a)At  [-Ksqs  +  Fs]  +  aAt[Ks+lqs+l  -  Fs+l  ]  =  0  (3.32) 

Using 

al  =  aAt  (3.33) 

a2  =  (l-a)At  (3.34) 

in  Eq.(3.32)  gives 

MAq-a2 [~Ksqs  +FS]  +  al[Ks+1qs+l  -Fs+1]  =  0  (3.35) 

Assuming  that  the  diffusivity  does  not  change  with  respect  to  time,  we  have 

[*„.]=[*,]  <3-36> 

Using  Eq.(3.36),  Eq.(3.35)  can  be  re-written  as 

MAq-a2 [-Ksqs  +  Fs]  +  al[Ksq,  +  KAq  - Fs+l ]  =  0  (3.37) 

Rearranging  to  bring  all  the  terms  involving  the  unknowns  to  the  left  side  gives 

[M  +  alKs  ]  A  q  =  -(al  +  a2)  [Ksqs  ]  +  a2[Fs]  +  al[Fs+1  ]  (3.38) 

Eq.(3.38)  is  solved  to  obtain  the  solution  for  the  s+1  time  step.  Therefore,  the  finite 
element  formulation  for  this  diffusion  model  can  be  described  by  the  following 
equations 

MAq  =  F  (3.39) 

where 

M  =[M  +a\Ks]  (3.40) 

F  =  -(al  +  a2)[Ksqs]  +  a2[Fs]  +  al[Fs+l]  (3.41) 
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3.3.7.  Boundary  conditions  for  diffusion  of  periodic  microstructures 

Periodic  configurations  can  be  analyzed  by  using  just  a  representative  volume  element 
(RVE)  or  unit  cell.  Similar  to  how  solid  mechanics  models  can  be  homogenized,  RVE 
models  of  periodic  microstructures  can  also  be  used  to  obtain  effective  diffusivities.  One 
noticeable  difference  with  the  solid  mechanics  models  described  in  the  earlier  sections  is 
that  they  deal  with  static  mechanics  whereas  the  diffusion  model  described  in  the 
previous  section  simulates  a  transient  behavior.  In  order  to  calculate  the  effective 
diffusivities,  the  concentration  distribution  in  the  model  at  steady-state  conditions  is 
required.  This  means  that  the  transient  part  of  Eq.(3.17)  is  omitted  making  it  a  static 
model. 

Once  the  unit  cell  is  chosen  for  the  periodic  microstructure,  certain  characteristics  can  be 
determined  based  on  the  fact  that  each  of  the  unit  cells  will  behave  in  the  same  manner. 
For  diffusion  at  steady-state,  the  periodic  conditions  state  that  the  concentration  gradient 
and  flux  distributions  are  identical  in  all  of  the  unit  cells.  This  can  be  expressed  as 


where  da  is  a  vector  of  periodicity  [43-44],  The  vector  of  periodicity  is  a  vector  from  a 
point  in  one  unit  cell  to  an  equivalent  point  in  an  adjacent  unit  cell. 

The  configurations  analyzed  in  this  work  are  in  general  heterogeneous  and  as  mentioned 
in  the  previous  section,  continuity  of  the  normalized  concentrations  is  imposed  in  order 
to  resolve  the  issue  of  discontinuous  concentrations  at  the  interface  of  two  different  base 
materials.  Therefore  all  the  formulations  and  models  are  defined  based  on  normalized 
concentrations,  C .  The  actual  concentrations  can  of  course  always  be  calculated  using 
Eq.(3.9).  In  some  ways  this  is  different  from  the  typical  homogenization  procedure  in 
solid  mechanics.  The  primary  variable  in  solid  mechanics  is  displacements  whereas  in 
diffusion,  the  typical  primary  variable  is  concentration,  which  is  generally  discontinuous 
across  different  base  materials.  This,  as  mentioned  earlier  necessitates  the  use  of 
normalized  concentrations,  which  is  continuous  across  different  base  materials.  To 
explain  the  subtle  differences  when  dealing  with  normalized  concentrations,  the 
procedure  to  determine  the  effective  diffusivity  properties  of  a  composite  with  circular 
fibers  in  a  periodic  square  array  is  described.  This  procedure  is  also  used  to  perform 
some  of  the  analyses  in  this  work. 
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The  approach  described  herein  is  consistent  with  Whitcomb  and  Tang’s  work[137]  but 
some  notations  have  been  changed  to  make  it  clearer.  Consider  a  discrete  unit  cell  of  a 
periodic  square  array  as  shown  in  Figure  3.3  and  assume  that  the  diffusing  mass  is 
macroscopically  flowing  in  the  horizontal  direction  and  therefore  there  is  no  flux  across 
the  top  and  bottom  edges.  Although  the  fiber  is  considered  to  be  impermeable  in  this 
work,  this  formulation  is  developed  assuming  that  both  the  matrix  and  fiber  are 
permeable  and  have  saturation  concentrations  of  Z)m,C"  and  Df,C J  respectively.  The 

matrix  is  assumed  to  be  isotropic  and  the  fibers  are  assumed  to  be  transversely  isotropic. 
Since  the  constituents  are  isotropic  in  the  transverse  plane  and  the  fibers  are  arranged  in 
a  square  array,  the  resulting  microstructure  will  have  the  same  effective  diffusivity  in  the 
x  and  y  directions,  denoted  by  Deff  .  Therefore,  in  order  to  obtain  the  effective 


Discrete  Unit 
Cell 


Equivalent 
Homogenized 
Unit  Cell 


Figure  3.3:  Boundary  conditions  for  the  discrete  unit  cell  and  the  equivalent 

homogeneous  unit  cell 


diffusivity  for  the  microstructure,  only  one  type  of  configuration  needs  to  be  analyzed 
with  an  imposed  concentration  gradient  in  the  x  direction.  Suppose  the  concentrations  on 
the  left  and  right  are  Clefi  and  Cnght,  respectively.  The  respective  normalized 
concentrations  are  obtained  by  dividing  the  concentrations  by  C" .  The  finite  element 
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model  of  the  configuration  can  be  analyzed  by  imposing  the  corresponding  normalized 
concentrations  on  the  left  and  right  edges.  The  results  will  show  a  continuous  variation 
of  the  normalized  concentration  across  the  domain  as  expected  but  the  actual 
concentrations  will  have  a  discontinuity  at  the  interface  between  the  fiber  and  the  matrix, 
if  they  have  different  saturation  concentrations.  It  is  convenient  to  define  an  effective 
property,  d  for  the  unit  cell  in  terms  of  volume  averaged  values  of  the  flux  in  the  x 

SC 

direction,  qx  and  the  normalized  concentration  gradient, - 

dx 

<3'45) 


where  the  angle  brackets  indicate  the  volume  average  of  the  bracketed  term. 

The  solution  can  be  post-processed  to  give  the  volume  averaged  flux  over  the  unit  cell. 
For  this  simple  geometry  and  boundary  conditions,  the  volume-averaged  normalized 
concentration  gradient  is  simply  (crtght  -  Cleft  ,  although  for  more  complicated 

models,  the  value  can  be  obtained  by  post-processing  the  solution. 

An  equivalent  homogenized  material  will  have  a  saturation  concentration  value  which  is 
the  volume-averaged  value  of  the  constituent  saturation  concentrations. 

c:=v„cz+vfc;  (3.46) 


In  the  corresponding  homogenized  unit  cell,  the  normalized  concentrations  at  the  right 
and  left  will  be  the  same  as  that  in  the  discrete  unit  cell  as  indicated  in  Figure  2.3.  The 
actual  concentrations  at  the  right  and  left  edge  in  the  homogenized  unit  cell  are  obtained 
by  using  Eq.(3.9).  Therefore  the  corresponding  concentration  on  the  left  and  right  will  be 

Q  left  right 

C;  — —  and  C;  — —  respectively  as  shown  in  Figure  3.3.  The  equivalent  concentration 

^ m 

gradient  can  be  written  as 


dC_ 

dx 
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(3.47) 


This  can  be  rewritten  in  terms  of  the  volume  averaged  normalized  concentration 

..  /dC\ 
gradients,  (  —  ) 


(3.48) 
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Combining  Eq.(3.45)  and  Eq.(3.48)  gives, 


which  gives  the  expression  for  effective  diffusivity  as 

Dff  =  — 

eJJ 

C h 


(3.49) 


(3.50) 


When  the  fiber  is  assumed  to  be  impermeable,  i.e.  Df  =  0 ,CJ  =  0 ,  the  expression  for 
the  effective  diffusivity  simplifies  to 


Deff 


vci 


(3.51) 


Under  such  an  assumption,  it  is  observed[12]  that  the  ratio - —  is  constant  for  a  fixed 

fiber  fraction,  regardless  of  the  value  of  the  matrix  diffusivity.  Let  this  ratio  be  defined 
by  the  following, 


D  = 


dc: 


(3.52) 


A  master  curve  can  be  obtained  showing  the  variation  of  D  with  fiber  fraction.  This 
master  curve  shown  in  Figure  3.4  is  valid  as  long  as  the  diffusion  follows  Fick’s  law. 
The  same  is  true  for  hexagonal  arrays  of  impermeable  fibers  and  Ref.  [12]  gives  a  simple 
curve  fit  for  both  master  curves.  This  makes  it  convenient  to  obtain  the  effective 
diffusivity  of  a  composite  with  impermeable  circular  fibers  for  various  fiber  fractions 
using  the  following, 


Deff 


DD 


V. 


(3.53) 


where  D  is  obtained  using  the  curve  in  Figure  2.4,  which  also  describes  the  formula  for 
the  curve  fit. 
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Figure  3.4:  Master  curve  for  impermeable  circular  fibers  in  square  array  showing 
variation  of  D  with  fiber  fraction  Vf  (refer  to  eq(3.52)  for  definition  of  D) 


3.4.  Oxidation 

This  section  describes  the  governing  equations  and  the  finite  element  formulation  for 
simulating  the  oxidation  behavior  in  pure  resin  and  in  composites.  The  oxidation  model 
can  be  considered  an  extension  of  the  diffusion  model  as  they  are  both  based  on  the 
same  conservation  law.  Some  of  the  issues  such  as  using  normalized  concentration  as  the 
primary  variable  are  common  to  this  analysis  as  well.  The  common  aspects  between  the 
two  analyses  will  be  highlighted  while  describing  the  special  circumstances  that  make 
this  analysis  different.  This  section  will  follow  the  same  format  as  the  section  on 
diffusion  and  starts  by  listing  the  equations  for  the  common  analysis  procedure.  This  is 
followed  by  the  derivation  of  the  weak  form  and  its  discretization  to  obtain  the  finite 
element  formulation. 

3.4.1.  Conservation  laws 

The  oxidation  process  in  a  polymer  is  a  combination  of  the  diffusion  of  oxygen  and  its 
consumption  by  reaction,  which  also  results  in  the  creation  of  by-products  such  as 
carbon  dioxide.  For  the  purposes  of  modeling  the  oxidation  of  polymers,  the  process  is 
assumed  to  be  dominated  by  the  diffusion  of  oxygen  into  the  polymer.  The  oxidation 
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model  that  is  used  in  this  research  effort  is  primarily  based  on  the  work  by  Pochiraju  et 
al[24-26]  in  which  they  used  the  conservation  of  mass  law  for  diffusion  with  a  term  to 
model  the  rate  of  consumption  of  the  diffusing  oxygen  during  chemical  reaction.  The 
governing  equation  can  be  expressed  as 


dc  dJi  D  n 

- + - !-  + j?  =  0 

dt  dx. 


(3.54) 


where  C  is  the  concentration  of  oxygen,  J ]  is  the  diffusion  flux  and  R  is  the  reaction 
rate  term. 

This  section  is  also  used  to  completely  define  the  reaction  rate  term  and  the  related 
quantities.  The  reaction  rate,  R  ?  in  general,  would  depend  on  the  concentration  of 
oxygen,  temperature  and  the  availability  of  un-oxidized  polymer.  As  the  oxygen  reacts 
with  the  polymer,  the  amount  of  polymer  available  for  oxidation  depletes  and  the  oxygen 
will  continue  to  diffuse  to  the  interior  of  the  polymer  to  react.  Depending  on  the  type  of 
polymer,  the  process  also  leads  to  a  reduction  in  the  molecular  weight  of  the  material 
due  to  chemical  bond  breakage  and  the  release  of  the  oxidation  by-products  [26].  The 
amount  of  polymer  available  for  oxidation  is  defined  by  an  oxidation  state  variable 

called  (f> .  The  value  of  the  oxidation  state  variable  at  which  the  polymer  is  considered  to 

be  completely  oxidized  with  no  more  polymer  available  for  reaction  is  defined  as  <j)ox. 

The  oxidation  state  can  be  physically  defined  to  be  the  ratio  of  the  current  weight  of  the 
material  over  its  original  un-oxidized  weight.  Therefore,  the  oxidation  state  $  has  a 

range  from  (j)ox  to  1  where  an  oxidation  state  value  of  1  denotes  the  un-oxidized  polymer. 

An  oxidation  state  value  between  (j)ox  and  1  indicates  that  the  material  is  partly  oxidized 

and  can  still  undergo  more  oxidation.  To  illustrate  this,  three  zones  were  defined  by 
Pochiraju  et  al[24-26]  as  shown  in  Figure  2.5.  Consider  that  the  left  end  of  the  idealized 
material  shown  in  the  figure  is  exposed  to  oxygen  and  the  oxidation  propagates  to  the 
right.  Zone  III  is  the  region  of  the  material  that  is  un-oxidized  with  an  oxidation  state  of 
1  and  as  the  oxidation  continues,  this  zone  gets  smaller  while  Zone  I  which  denotes  the 
fully  oxidized  material  with  an  oxidation  state  of  (j)ox  increases.  The  zone  in  between 

where  the  oxidation  state  is  between  (j)ox  and  1  is  called  the  active  zone  and  is  denoted 

by  Zone  II.  The  expression  for  calculating  the  oxidation  state  variable  is  described  later 
in  this  section. 
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Figure  3.5:  Oxidation  zones  and  corresponding  values  of  the  oxidation  state  variable 


When  (j>  =  <j)ox  at  a  material  point,  R  =  0  and  the  process  simplifies  to  just  diffusion  at  that 

point.  It  is  assumed  that  the  effects  of  (j>,T  and  C  on  R  are  separable  such  that  the 
function  R  can  be  expressed  as 


R  = 


/  (C)R0(T) 


(3.55) 


R0{T )  is  the  “saturated  reaction  rate”  (as  described  in  Ref. [26]),  which  describes  the 

dependence  of  the  reaction  rate  on  the  temperature  under  saturation  conditions.  The 
polymer  is  considered  saturated  when  it  has  the  maximum  amount  of  oxygen  possible 
for  the  given  temperature  and  pressure.  The  leading  factor  and  /(C)  in  the  expression 
both  have  a  range  from  0  to  1 .  The  leading  factor  models  the  dependence  of  the  reaction 
rate  on  the  availability  of  polymer  that  can  be  oxidized  such  that  R  is  maximum  when  $ 

has  a  value  of  1  and  linearly  decreases  to  zero  when  (f>  =  (j>gx .  The  function  /(C)  models 

the  dependence  of  the  reaction  rate  on  the  oxygen  concentration.  For  modeling  oxidation 
in  polyimide  resin  systems  like  PMR-15  as  implemented  by  Pochiraju,  the  function  /(C) 
is  taken  from  the  work  by  Colin  et  al[46-47]. 


/(C)  = 


2pC_ 

1  +pc 


1-  ^ 

2(1  +  pC) 


(3.56) 


The  value  of  (3  is  determined  by  using  weight  loss  data  obtained  from  specimens  aged 
at  two  different  oxygen  partial  pressures  i.e.  at  two  different  saturation  conditions, 
typically  in  pure  oxygen  and  air.  The  details  of  this  procedure  are  given  in  Ref.  [26],  The 
following  ratio  is  obtained  from  the  experimental  work  by  Abdeljaoued[40], 


weight  lossair 
weight  loss 

o  pure  oxygen 


(3.57) 
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Assuming  that  the  weight  loss  is  proportional  to  the  reaction  rates,  the  ratio  of  the  weight 
loss  from  the  two  specimens  would  be  the  same  as  the  ratio  of  the  reaction  rates  and 
would  give  the  following  equation, 


/?(C  =  0.79 mo// m3 , 288°C)  _  /^(288°C)/(C  =  0.79mo//m3) 
R(C  =  3.74mo//m3 ,288°C)  ~  /^(288°C)/(C  =  3.74 mo//m3) 


This  can  be  solved  to  obtain  three  roots  for  (3  of  which  only  one  is  non- zero  and  has  a 
value  of  0.919. 


For  modeling  neat  PMR-15  resin,  <f>m  has  a  value  of  0.187,  which  is  taken  from 

Pochiraju’s  work  [24-26],  This  value  is  determined  from  experimental  weight  loss  data 
and  the  method  is  described  in  Ref.  [26].  The  oxidation  state  variable  can  be  related  to 
the  weight  loss  of  the  material  as  follows 


dd>  dW 

— —  oc - 

dt  dt 


(3.59) 


where  W  is  the  weight  of  the  material. 

Assuming  that  the  rate  of  change  of  weight  is  proportional  to  the  reaction  rate  gives, 


dW 

dt 


cc  —R 


Combining  Eq.(3.59)  and  Eq.(3.60)  gives  the  following. 


d(j) 

dt 


=  -aR 


(3.60) 


(3.61) 


where  a  is  a  proportionality  parameter  that  is,  in  general,  time  and  temperature 
dependent. 

Using  Eq.(3.61),  the  following  expression  for  calculating  (j)  can  be  obtained 


tj)  =  max 


t 

l-\a(g)R(g)dg 

0 


(3.62) 


An  issue  that  arises  when  analyzing  oxidation  in  heterogeneous  materials  is  that 
although  the  oxidation  state  value  for  any  material  has  an  upper  limit  of  1,  its  lower  limit 
for  different  materials  is  not  necessarily  the  same.  This  makes  it  inconvenient  to  make 
comparisons  as  to  how  much  oxidation  has  taken  place.  For  example,  the  same  oxidation 
state  value  for  two  different  materials  need  not  imply  that  they  are  equally  close  to  being 
fully  oxidized  or  that  they  have  the  same  amount  of  material  left  to  oxidize.  In  order  to 
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make  this  comparison  easier,  a  new  variable  is  introduced  called  the  oxidation  level 
denoted  by  O  .  The  oxidation  level  variable  linearly  scales  the  oxidation  state  variable 
so  that  all  materials  have  an  oxidation  level  that  ranges  from  0  to  1.  This  relation  is  given 
by, 


o  = 


(3.63) 


For  the  same  reasons  described  in  Section  3.3.1,  the  differential  equation  described  in 
Eq.(3.54)  is  re-written  in  terms  of  normalized  concentrations, 


8t  8x: 


(3.64) 


3.4.2.  Concentration  gradients 

Just  as  in  the  diffusion  analysis,  the  oxidation  analysis  uses  simple  gradients  of  the 


8C 

concentration,  — 

8xi 


3.4.3.  Constitutive  relations 

The  relationship  between  flux  and  the  concentration  gradient  is  given  by  Fick’s  first  law, 


J:  =  -D, 


8C_ 

8x: 


(3.65) 


where  Z>  is  the  2nd  order  diffusivity  tensor.  The  Latin  subscripts  i  and  j  denote  the 

coordinate  direction  and  range  from  1  to  3  for  a  three  dimensional  formulation.  The 
constitutive  relationship  is  different  from  that  in  the  diffusion  analysis  described  in 
Section  3.3.3  and  that  is  because  the  diffusivities  of  the  un-oxidized  and  oxidized 
polymer,  in  general,  will  be  different.  The  diffusivity  of  the  polymer  in  the  active 
oxidizing  zone  (where  (j)ox  <  (j)  <  1 )  is  assumed  to  vary  linearly  between  the  un-oxidized 

polymer  diffusivity  and  the  fully  oxidized  polymer  diffusivity  and  is  given  the  following 
expression 


D«=(D:), 


1-d 


+ 


A)„ 


OX  J 


\z± 

Wc 


(3.66) 


ox  J 


Again,  Eq.(3.65)  is  re-written  in  terms  of  the  normalized  concentration 
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3.4.4. 


Boundary  conditions 


The  boundary  conditions  are  defined  in  the  same  manner  as  the  diffusion  analysis.  The 
flux  boundary  conditions  are  given  by 

q  =  -niJi  on  S  (3.68) 

And  the  normalized  concentration  boundary  conditions  are  given  by 

C  =  C  on  5  (3.69) 

where  C  is  the  specified  normalized  concentration  on  the  boundary  5 


3.4.5.  Weak  form 


The  same  procedure  is  applied  as  described  in  Section  3.3.5  for  the  derivation  of  weak 
form  for  diffusion  analysis.  The  difference  is  in  the  inclusion  of  the  reaction  rate  term  in 
the  case  of  the  oxidation  analysis.  The  equation  of  conservation  mass  is  multiplied  by  a 
variation  of  the  normalized  concentration  and  integrated  over  the  volume  to  obtain  the 
weighted  integral  form, 


—  +  —  J.+R  \dV=0 

dt  dx: 


(3.70) 


where  SC  is  an  arbitrary  variation  of  the  normalized  concentration. 
Integration  by  parts  gives  the  weak  from, 


Using  Eq.(3.67) 
formulation, 


S'CC'—-C-^-Ji+RSC 

dt  dx. 


dV  +  jSCniJidS=  0  (3.71) 


and  Eq.(3.68)  in  Eq.(3.71)  gives  the  basis  for  the  finite  element 


6C  dSC 

- 1 - 

dt  dxt 


8C 

CD..  — —  +  RSC 

dx: 


\dV  = 


s 


(3.72) 


3.4.6.  Discretization  of  weak  form  and  time  integration 

Again,  the  same  basic  procedure  is  applied  as  described  in  Section  3.3.6  for  the 
derivation  of  finite  element  formulation.  On  the  other  hand,  there  are  some  details  that 
are  quite  different  from  the  diffusion  analysis.  This  is  because  of  the  reaction  rate  term 
and  the  non-linear  expression  of  the  diffusivity  in  the  weak  form. 

Over  a  typical  finite  element,  the  normalized  concentration  is  approximated  by 

C(x,t)  =  Na(x)Ca(t)  (3.73) 


where  Na  are  the  interpolation  functions  and  Ca  are  the  nodal  normalized 

concentrations.  The  subscripts  with  Greek  letters  range  from  1  to  the  number  of 
interpolation  functions. 

After  discretizing  the  weak  form  using  Eq.(3.73)  and  SC  =  NaSCa ,  the  following  set  of 
equations  are  obtained, 


dC„  dN 
N  rwN  _ £  + 

8  '  dt 


dN „  - 

— —  C°D..  — —  CB  +  NR 

dx  ,J  dx :  P  a 


dv  =  J  NjdS 


In  matrix  form  this  can  be  written  as 

M  a/A  ft  +  KapClp  +  Ra  =  Fa 


where 


M 


J[iV„C"iV„]rfV 


dL 

dx. 


CmD 


w 


K=\[n„RYv 

V 

F«=\NJdS 


(3.74) 

(3.75) 

(3.76) 

(3.77) 

(3.78) 

(3.79) 
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qa=Ca  and  qa=^% 


(3.80) 


Just  as  in  the  case  of  the  diffusion  analysis,  an  approximation  is  used  for  the  time 
derivative  term  in  Eq.(3.74)  whereby  the  solution  at  a  particular  instant  in  time  is 
determined  based  on  the  solution  history.  The  same  ‘time-marching’  procedure  from 
Section  3.3.6  is  used  to  numerically  solve  Eq.(3.74).  Using  Eq.(3.75),  the  following 
equations  can  be  written  for  two  consecutive  time  steps,  t=ts  and  t=ts+i, 

Mqs  +  Ksqs  +R*  —Fs  =  0  (3.81) 

Mqs+l  +  Ks+1qs+1  +  Rs+1  -  Fs+1  =  0  (3.82) 

Using  the  alpha  family  of  approximations  [45]  gives 

(1  -  a)qs  +  aqs+l  =Aq/  At  (3.83) 

Multiplying  Eq.(3.83)  by  AtM  gives 


(1  -  a)AtMqs  +  aAtMqs+1  =  MAq  (3.84) 

Rearranging  the  terms  in  Eq.(3.84)  gives  an  expression  for  aAtMqs+l 

aAtMqs+l  =  MAq  -  (1  -  a)AtMqs  (3.85) 

Multiplying  Eq.(3.82)  throughout  by  aAt  gives 

aAtMqs+l  +  aAt  [_Ks+iqs+l  +  R*s+l  -  Fs+l  ]  =  0  (3.86) 

Substituting  Eq.(3.85)  into  Eq.(3.86)  gives, 

MAq  -  (1  -  a)AtMqs  +  aAt  [_Ks+xqs+x  +  R*s+1  -  Fs+1  ]  =  0  (3.87) 

And  substituting  for  Mqs  from  Eq.(3.81)  in  Eq.(3.87)  gives  the  governing  equations 

M^-(\-a)M[-KA-R]+F,}+aM  +  R'„,  -  FM ]  =  0  (3.88) 

Using  Eq.(3.33)  and  Eq.(3.34)  in  Eq.(3.88)  gives 

MA9-«2[-;:A-fi;+F,]+al[*:„,9„,+«;,-F„,]  =  0  (3.89) 
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A  Taylor  Series  expansion  is  used  on  the  terms  in  Eq.(3.89)  that  are  dependent  on  the 
unknown  solution,  ( qs+l ),  which  are  Ks+iqs+1  and  R*s+l .  Indices  will  be  used  in  the  next 

few  steps  in  order  to  make  the  operations  involved  clear.  Ignoring  the  higher  order  terms 
in  the  Taylor  Series  expansion  gives  the  following  expression, 


(«/>)„, +KL  +K), +: 


d 


MM 


dqc 


d{K) 

Aq  +^—^Aq 

dqg 


(3.90) 


The  partial  derivatives  in  the  expression  above  are  very  complex  and  therefore  the  aim  is 
to  obtain  an  approximation  for  the  expression.  It  is  assumed  that  for  sufficiently  small 
time  steps,  the  error  is  minimal  and  certain  approximations  can  be  made.  Similar 
approximations  have  been  made  in  Pochiraju’s  oxidation  model  [24-26],  One  approach 

to  obtain  an  approximate  expression  for  is  by  assuming  that  only  /(C)  from 

8qp  ' 

Eq.(3.55)  depends  on  C .  This  would  make  it  simpler  to  take  a  derivative  of  R*a  with 


respect  to  the  nodal  variables,  qg .  This  approach  will  be  evaluated  in  future  work  but  for 
this  work,  it  is  assumed  that  if  the  time  step  is  sufficiently  small  that  (// )  ^  in 

Eq.(3.90)  can  be  approximated  by  (R* )  (or  mathematically,  ^(^1  =  o).  The  remaining 

v  7' 


partial  derivative  in  Eq.(3.90)  can  be  expressed  as 


dqc 


dqc 


M, 


dqc 


(3.91) 


The  term  ^  is  not  convenient  to  compute  because  Kap  depends  on  /  ,  which  is  a 

complex  function  of  the  solution  (see  Eq.(3.62)).  Again,  it  is  assumed  that  for 

d(Ka0 ) 

sufficiently  small  time  steps,  Eq.(39)  can  be  approximated  by  assuming  — — —  =  0. 

dq( 

Thus  Eq.(3.91)  becomes 


dq. 


(3.92) 
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Therefore  Eq.(3.90)  becomes 

MJ«*L+teL («A+Kl <3-93> 

Substituting  Eq.(3.93)  in  Eq.(3.89)  gives 

M Ag -  a.2 [~K,q,  -  R]  +  F, ]  +  al [( K,q,  +  R])  +  K, Ag - F„, ]  »  0  (3,94) 

Rearranging  to  bring  all  the  terms  involving  the  unknowns  to  the  left  side  gives 

[M  +  alKs  ] A q  *  -(al  +  al) [Ksqs  +  R*s  ]  +  a2 [Fs] +  al [F+l  ]  (3.95) 

Eq.(3.95)  is  solved  to  obtain  the  solution  for  the  s+1  time  step.  For  sufficiently  small 
time  steps,  it  is  seen  that  this  approximation  is  reasonable  because  a  parametric  study 
with  different  time  step  sizes  showed  the  model  appearing  to  converge  to  the  same 
solution.  Therefore,  the  finite  element  formulation  for  this  oxidation  model  can  be 
described  by  the  following  equations 

MAq  =  F  (3.96) 

where 

M  =[M  +a\Ks]  (3.97) 

F  =  -(al  +  a2)[Ksqs  +  R*s  ]  +  a2  [Fs  ]  +  al  [Fs+l  ]  (3.98) 


To  arrive  at  this  formulation,  a  number  of  approximations  were  made  to  simplify  the 
nonlinearity.  Typically,  when  solving  a  nonlinear  equation  numerically,  a  ‘residual’  is 
driven  to  zero  by  iterating.  In  this  implementation,  there  is  no  iterating  at  each  time  step 
in  order  to  drive  a  ‘residual’  to  zero.  This  is  because  it  was  found  that  the  even  without 
iterating,  the  results  were  found  to  be  reasonably  close  to  that  from  Pochiraju’s  model. 

An  important  part  of  the  oxidation  analysis  is  post-processing  the  results  of  the 
simulation  to  provide  a  measure  of  the  oxidation  behavior.  The  oxidation  behavior  is 
visualized  in  terms  of  the  growth  of  the  oxidation  layer.  The  oxidation  layer  initiates 
from  the  surfaces  exposed  to  the  oxygen  and  grows  into  the  interior  as  the  material 
becomes  oxidized.  Although  ideally  the  material  is  said  to  have  started  oxidizing  when 
the  oxidation  level  drops  below  1,  the  oxidation  layer  thickness  is  defined  by  the  point  at 
which  the  oxidation  level,  O,  dips  below  0.99,  indicating  that  1%  of  the  oxidizable 
material  has  oxidized.  Therefore,  an  element  is  assumed  to  have  started  oxidizing  if  the 
oxidation  level  at  each  of  the  material  integration  points  falls  below  0.99.  If  the 
oxidation  state  is  above  0.99,  the  element  is  assumed  to  be  un-oxidized  and  if  it  is  below 
0.01  it  is  assumed  to  be  fully  oxidized.  A  post-processing  routine  was  written  that 
calculated  the  growth  of  the  oxidation  layer  in  the  ID  model.  This  involved 
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extrapolating  the  oxidation  state  values  from  the  integration  points  to  the  nodal  points, 
averaging  the  extrapolated  values  at  a  node  if  the  node  shared  elements  of  the  same 
material  and  solving  for  the  location  in  the  model  where  the  oxidation  level  value  met 
the  specified  upper  and  lower  limits.  This  routine  was  also  generalized  to  work  for  2-D 
and  3-D  models.  Note  that  the  prescribed  upper  and  lower  limits  of  0.99  and  0.01 
respectively  are  valid  only  for  a  completely  oxidizable  material  such  as  neat  resin.  When 
dealing  with  homogenized  material  such  as  a  tow,  the  entire  material  does  not  oxidize 
because  the  fibers  are  assumed  to  be  inert  and  therefore  the  prescribed  limits  will  be 
different.  In  such  a  case,  the  upper  limit  that  defines  the  oxidation  layer  thickness  is  the 
oxidation  level  when  1%  of  the  resin  in  a  homogenized  tow  material  point  is  oxidized. 
This  upper  limit  is  given  by 

O  =  1-0.01V  (3.99) 

The  model  described  in  this  section  was  implemented  and  the  oxidation  layer  growth 
was  simulated  for  neat  PMR-15  resin  using  a  ID  model.  The  results  were  compared  with 
that  from  Pochiraju’s  simulation.  For  a  200-hr  simulation,  it  was  found  that  both  models 
agree  closely  in  predicting  the  Zone  I  thickness.  The  difference  is  negligible  in  the 
beginning  of  the  simulation  and  grows  to  a  maximum  difference  at  200  hours  when  the 
current  model  predicts  a  thickness  107  microns  compared  to  104  microns  predicted  by 
Pochiraju’s  model.  Both  models  predict  the  Zone  II  thickness  to  be  almost  constant 
throughout  the  200  hr  simulation.  Pochiraju’s  model  gives  a  Zone  II  thickness  of  19 
microns  while  the  current  model  under  predicts  by  21%  with  a  value  of  15  microns.  The 
cause  of  this  difference  could  be  the  implementation  of  the  two  models.  Pochiraju’s  ID 
model[26]  uses  a  modified  implementation  of  odel5s  and  Pdepe  solvers  in  MATLAB  to 
solve  the  governing  equation,  Eq.(3.54).  The  current  model  on  the  other  hand  uses  a 
standard  one-dimensional  finite  element  implementation  based  on  Eq.(3.96).  For  the 
purposes  of  investigating  the  effect  of  oxidation  on  the  mechanical  response  of  the 
composites  using  this  material  system,  it  is  assumed  that  the  thickness  of  Zone  I  alone  or 
the  overall  thickness  (Zone  I  +  II)  that  is  of  primary  concern.  Thus,  if  the  overall 
thickness  is  considered,  the  difference  between  the  two  models  is  around  21%  in  the 
beginning  and  drops  to  about  2%  at  200  hours,  which  is  assumed  to  be  negligible  for  the 
purposes  of  this  particular  research  effort. 

The  various  material  input  properties  required  for  specifying  the  equations  in  the 
oxidation  model  are: 

1.  The  diffusivities  for  the  oxidized  and  un-oxidized  material,  Dox ,  Dunox 

2.  Saturated  reaction  rate,  R0 

3.  Dependence  of  reaction  rate  on  concentration,  /(C)  and  the  constant  f3 

4.  Value  of  oxidation  state  when  fully  oxidized,  (/>ox 

5.  Weight-reaction  proportionality  parameter,  a 
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3.4.7.  Boundary  conditions  for  oxidation  of  periodic  microstructures  and 

homogenization 

The  oxidation  response  in  polymers  and  PMCs  is  a  nonlinear  transient  behavior.  Just  as 
in  the  case  of  diffusion,  in  order  to  simulate  oxidation  for  periodic  microstructures,  the 
transient  part  of  the  behavior  needs  to  be  removed  effectively  looking  at  the 
microstructure  at  steady-state  conditions.  Under  oxidation  behavior,  steady-state 
conditions  imply  that  all  the  material  in  the  microstructure  is  oxidized.  But  when  all  the 
material  is  oxidized,  the  behavior  reverts  to  typical  steady-state  diffusion  behavior. 
Therefore,  it  is  not  practical  to  obtain  the  effective  oxidation  material  properties  in  this 
manner.  Instead,  other  strategies  are  explored  for  simulating  oxidation  in  periodic 
microstructures. 

In  order  to  model  oxidation  in  a  textile  composite,  it  is  necessary  to  obtain  effective 
properties  for  the  tows  because  it  is  impractical  or  even  impossible  to  discretely  model 
the  entire  microstructure.  This  section  will  describe  approaches  for  obtaining  effective 
oxidation  material  properties  for  tows. 

When  trying  to  replace  a  heterogeneous  material  with  a  homogenized  material  in  model, 
it  is  reasonable  to  assume  that  some  or  possibly  all  of  these  properties  might  change.  It 
can  also  be  expected  that  there  could  be  more  than  one  way  that  these  properties  can 
change  to  demonstrate  the  same  overall  behavior  as  a  discretely  modeled  heterogeneous 
microstructure.  There  are  at  least  two  approaches  for  achieving  this  goal.  One  is  to  use  a 
multi-scale  analysis  that  keeps  track  of  the  ‘local’  information  such  as  oxidation  state 
and  actual  average  concentration  in  the  constituent  matrix  in  the  homogenized  material. 
Given  this  information,  it  would  be  possible  to  calculate  the  reaction  rate  R  at  a 
particular  material  point  in  the  tow’s  constituent  matrix  using  Eq.(3.55).  Next,  the 
effective  reaction  rate  for  the  larger  scale  homogenized  tow  is  determined  by  a  simple 
rule  of  mixtures  and  plugged  into  the  governing  equations.  When  the  equations  for  a 
time  step  are  solved,  the  calculated  concentrations  are  transformed  back  to  the  local 
scale  using  a  rule  of  mixtures  in  order  to  keep  track  of  the  oxidation  state  of  the 
constituent  matrix.  Thus,  a  continuous  transfer  of  information  between  the  two  scales 
needs  to  be  maintained  throughout  the  simulation.  For  this  work,  another  approach  is 
used  where  effective  oxidation  properties  for  the  homogenized  material  are  determined 
thereby  eliminating  the  need  to  go  back  and  forth  between  the  two  scales.  A  few 
assumptions  are  made  in  order  to  determine  the  effective  material  properties,  Dox,Dunox, 

r0  ,  /(C) ,  fi  ,<f>ox  and  a.  These  assumptions  and  the  procedure  to  determine  the  properties 
are  described  in  the  remainder  of  this  section. 

In  this  work,  the  fibers  in  the  tows  are  idealized  to  be  in  a  square  array  and  the  fibers  are 
assumed  to  be  impermeable  and  do  not  oxidize.  While  there  are  other  factors  that  can 
influence  the  oxidation  behavior  in  composites  such  as  the  properties  of  the  fiber/matrix 
interface  or  interphase,  they  are  not  taken  into  account  for  the  homogenization  model 
described  in  this  work.  Cracks  in  the  matrix  or  along  the  fiber/matrix  interface  can  also 
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affect  the  oxidation  behavior  by  allowing  rapid  ingress.  Depending  on  the  type  of 
damage  that  is  observed  in  these  composites,  it  might  be  possible  to  account  for  their 
effects  in  the  homogenized  model.  For  example,  if  the  damage  is  diffuse,  the 
homogenized  properties  can  be  degraded  appropriately  or  if  the  damage  is  confined  to 
certain  areas,  cracks  can  be  introduced  in  the  homogenized  model.  Since  this  model  does 
not  consider  factors  such  as  the  effect  of  damage,  fiber/matrix  interface  or  cracks  in  the 
composites,  the  only  oxidation  material  property  that  is  different  when  considering  axial 
and  transverse  growth  is  the  diffusivity.  The  axial  diffusivity  is  largely  governed  by  a 
rule  of  mixtures  and  exhibits  simple  behavior.  Therefore,  particular  attention  is  not  paid 
to  the  axial  oxidation  behavior.  Moreover,  in  realistic  applications,  the  surfaces  exposed 
to  oxidizing  environments  are  mostly  parallel  to  the  fibers.  The  laminate  configurations 
that  are  analyzed  in  this  work  are  chosen  based  on  these  considerations  and  therefore,  the 
oxidation  behavior  depends  on  the  transverse  oxidation  material  properties. 


3.4.8.  Oxidation  parameters 

This  sub  section  will  describe  some  of  the  effective  material  properties  mentioned  in  the 
previous  sub  section  related  to  boundary  conditions  for  oxidation  of  periodic 
microstructures. 

3.4.8.1.  Diffusivities  for  the  oxidized  and  un-oxidized  material  Dox,Dunox 

The  diffusivities  on  its  own  only  define  the  mass  flow  of  oxygen  in  the  material.  It  will 
be  assumed  that  the  oxidation  state  continues  to  have  a  linear  effect  on  the  effective 
diffusivities  of  the  homogenized  material.  The  effective  diffusivity  can  be  determined  by 
just  modeling  the  diffusion  without  the  need  for  modeling  the  oxidation  behavior.  The 
procedure  for  determining  effective  diffusivity  as  described  in  Section  3.3.7  is  used  to 
obtain  the  effective  diffusivities  for  the  oxidized  and  un-oxidized  material. 

3.4.8.2.  Saturated  reaction  rate,  R0 

Since  the  matrix  is  the  only  material  that  is  oxidizing,  the  effective  saturated  reaction 
rate  would  be  expected  to  be  related  to  the  amount  of  matrix  in  the  unit  cell.  It  is 
assumed  that  the  relationship  follows  a  rule  of  mixtures  (with  the  fiber  having  a  reaction 
rate  of  zero).  That  is 

^  =  V„K"nx  (3-100) 

3.4.8.3.  Dependence  of  concentration  on  reaction  rate  on  f(C )  and  (3 

The  term  /(C)  models  the  dependence  of  the  reaction  rate  on  the  oxygen  concentration. 
Colin’s  expression[46-47]  given  in  equation  (3.56),  which  is  used  as  /(C)  to  model  the 
neat  PMR-15  polymer  will  be  used  for  the  homogenized  tow  as  well.  It  is  assumed  that 
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the  same  expression  models  the  behavior  of  the  homogenized  tow.  As  mentioned  in  the 
previous  section,  to  determine  the  value  of  (5  the  ratio  of  the  weight  loss  of  the  material 
at  two  different  saturation  conditions  is  needed.  Due  to  lack  of  access  to  experimental 
data  on  oxidation  weight  loss  for  unidirectional  laminates,  it  is  assumed  that  the  ratio  of 
weight  loss  is  the  same  as  that  of  the  neat  PMR-15  polymer.  Therefore,  /?  has  the  same 
value  as  that  of  the  neat  resin,  which  is  0.919. 

3.4.8.4.  Oxidation  state  when  fully  oxidized,  (j)ox 

This  is  a  property  that  is  calculated  by  determining  the  weight  loss  of  the  material  when 
oxidized.  For  example,  a  value  of  0.2  implies  that  the  fully  oxidized  material  weighs 
about  20%  of  its  original  weight.  Due  to  lack  of  access  to  weight  loss  data  for 
unidirectional  laminates,  an  estimate  for  (j)ox  is  obtained  based  on  the  assumption  that  the 

fiber  does  not  lose  weight  during  oxidation.  Therefore,  the  effective  (f>ox  would  be  given 
by 

(3.101) 


3.4.8.5.  Weight-reaction  proportionality  parameter,  a 

In  general,  the  proportionality  parameter  a  is  time  and  temperature  dependent.  The 
value  of  a  for  the  neat  resin  is  determined  by  examining  the  oxidation  layer  growth.  On 
comparison  of  the  simulation  results  with  the  experimental  results,  Pochiraju  [9]  found 
that  the  oxidation  behavior  was  better  simulated  when  the  proportionality  parameter  was 
linearly  decreased  over  time  from  0.01  to  0.0033  for  the  first  40  hours  of  oxidation  and 
then  remains  constant  at  0.0033.  The  value  of  a  for  the  homogenized  tow  is  assumed  to 
follow  the  same  as  that  of  the  neat  PMR-15  resin. 

3.5.  Coupled  mechanical-oxidation  analysis 

A  coupled  mechanical-oxidation  analysis  model  was  developed  to  predict  damage 
initiation  and  progression  in  textile  composites  under  an  oxidizing  environment. 
Although  the  analyses  performed  in  this  work  assumes  only  one-way  coupling,  the 
underlying  analysis  model  forces  no  such  restriction  and  can  account  for  full  coupling 
between  the  mechanical  and  oxidation  analysis.  This  section  describes  the  coupled 
analysis  model  used  in  this  work  followed  by  the  constitutive  relations  used  to  the 
couple  the  two  analyses. 

One  component  of  the  coupled  analysis  is  the  oxidation  analysis  that  simulates  the 
diffusion  of  oxygen  into  the  composite  and  tracks  how  much  the  material  has  oxidized. 
The  second  component  is  the  damage  progression  analysis  that  can  track  the  damage  in 
the  material  and  degrade  the  properties  of  the  damaged  regions.  The  theory  and  finite 
element  formulation  behind  both  the  oxidation  analysis  and  the  damage  progression 
analysis  is  provided  in  the  previous  sections  and  they  are  adapted  to  use  in  this  coupled 
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analysis  model.  The  coupling  between  the  two  analyses  is  enabled  by  constitutive 
relations.  The  full  coupling  requires  a  constitutive  relation  relating  the  oxidation  state  to 
the  mechanical  properties  and  another  constitutive  relation  relating  the  mechanical  state 
to  the  oxidation  material  properties.  In  this  work,  all  the  configurations  that  were 
analyzed  assumed  only  a  one-way  coupling  with  the  oxidation  state  affecting  the 
mechanical  properties  of  the  model.  The  remainder  of  this  section  describes  the 
algorithm  for  this  one-way  coupled  model.  The  implementation  and  the  details  involved 
with  performing  an  actual  coupled  analysis  on  a  textile  composite  configuration  will  be 
explained  in  a  later  section. 

Since  the  analysis  assumes  only  one-way  coupling  and  the  mechanical  state  does  not 
affect  the  oxidation  material  properties,  the  oxidation  analysis  can  be  performed 
independent  of  the  damage  progression  analysis.  Therefore,  the  output  from  the 
oxidation  simulation  can  be  used  by  the  model  as  needed.  The  oxidation  analysis  output 
contains  the  oxidation  state  distribution  in  the  configuration  at  different  times  in  the 
simulated  oxidation  time  period.  The  damage  progression  analysis  described  is  a  quasi¬ 
static  analysis  where  the  loading  is  ramped  up  but  is  essentially  time-independent.  On 
the  other  hand,  all  the  coupled  models  analyzed  in  this  work  assume  a  constant 
mechanical  loading  while  the  configuration  is  undergoing  oxidation.  Therefore,  the 
damage  progression  analysis  cycles  through  each  of  the  time  data-points  in  the  simulated 
oxidation  time  period  and  performs  the  following  steps  -  Load  the  oxidation  state  for  the 
particular  time  data-point,  modify  the  mechanical  properties  and  iterate  to  converge 
upon  the  final  damage  state  for  the  corresponding  time  data-point.  This  is  illustrated  in 
the  flowchart  for  the  algorithm  shown  in  Figure  3.6. 

3.6.  Constitutive  Relations 

Experimental  results  show  that  oxidation  causes  damage  in  the  oxidized  material  which 
can  ultimately  affect  the  mechanical  properties  of  the  composite  [26].  Oxidation  is  found 
to  affect  the  mechanical  properties  of  fibers  [26].  But  it  is  not  trivial  to  characterize  the 
damage  and  its  effects  on  the  mechanical  properties  of  the  composites.  The  underlying 
mechanisms  and  the  properties  of  the  fiber/matrix  interface  and  interphase  have  not  been 
fully  understood  yet.  Shrinkage  of  the  matrix  due  to  oxidation  is  theorized  to  be  among 
the  factors  causing  delaminations  on  the  fiber  matrix  interface  [30].  These  cracks  can 
further  affect  the  oxidation  behavior  by  allowing  oxygen  to  penetrate  the  material  faster. 
But  the  effects  of  the  mechanical  or  physical  damage  on  the  oxidation  behavior  are  not 
being  considered  in  the  simulations  used  in  this  work.  This  section  will  describe  the  type 
of  constitutive  relations  used  in  the  simulations  that  were  performed  in  this  work. 

This  constitutive  relation  or  degradation  scheme  is  similar  in  some  respect  to  the 
property  degradation  scheme  based  on  mechanical  damage.  They  are  similar  in  the  sense 
that  the  engineering  moduli  are  modified  to  account  for  the  effect  of  the  oxidation.  The 
constitutive  relation  quantifies  the  amount  of  damage  in  terms  of  strength  and  stiffness 
degradation  based  on  the  oxidation  level  of  the  material  in  the  composite  (see  Eq.(3.63)). 
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Let  EtJ ,  Gtj  and  vn  be  the  original  extensional  moduli,  shear  moduli  and  Poisson’s 

ratio  respectively  and  Ey ,  Gy  and  Oy  be  the  degraded  extensional  moduli,  shear 

moduli  and  Poisson’s  ratio  respectively.  6  (i  =  1..9)  are  the  percentages  by  which  the 

nine  engineering  properties  change  when  the  material  is  completely  oxidized.  Remember 
that  the  oxidation  level,  <J>  ranges  from  1,  which  means  un-oxidized  to  0,  which  means 
fully  oxidized.  The  amount  of  degradation  is  assumed  to  vary  linearly  with  the  oxidation 
level.  Therefore,  a  typical  property  degradation  scheme  will  look  like: 

^  =  (l  +  (l_0)Wl,  ~En  =  (l  +  (l-<t>)b2)E22,  £^  =  (l  +  (l-0)&3)£33 

G^  =  (l  +  (l-<D)h4)G12,  G~  =  (1  +  (1-0)65)G23,  G^  =  (1  +  (1  - 0)h6)G33  (3.102) 

l>i2  —  (1  l  (1  —  0)67  )r>12 ,  u23  —  (1  +  (1  —  0)6g  )u23 ,  u13  —  (1  +  (1 —  ®)b^  )t^33 

For  example,  if  En  is  reduced  by  20%  when  the  material  is  fully  oxidized,  then 
bx  =  -0.2 .  If  the  En  property  needs  to  be  degraded  for  a  material  that  has  an  oxidation 
level  of  0.3,  the  new  modulus,  according  to  Eq.(3.102)  would  be  given  by 

=  (1  +  (1  —  0)6,  )En  =  (1  -  0.7  x  0.2)£u  =  0.86^!  (3. 103) 

Note  that  in  this  general  framework,  the  diagonal  as  well  as  non-diagonal  entries  of  the 
compliance  matrix  can  be  affected  independently. 

Similarly,  the  strength  can  also  be  degraded  based  on  the  amount  of  oxidation  the 
material  has  undergone.  In  this  work,  the  strengths  under  compression  are  assumed  to  be 
the  same  as  the  strengths  under  tension.  Let  Sn  (i  =  1..6)  denote  the  original  strengths  of 

the  material  in  the  different  stress  components  (in  Voigt  notation)  and  St,  (i  =  1..6)  be 
the  degraded  strengths.  Let  the  strength  degradation  parameters,  dt,(i- 1..6)  be  the 

corresponding  factors  by  which  the  strengths  would  be  degraded  if  the  material  was  fully 
oxidized.  Again,  a  linear  dependence  on  the  oxidation  level,  <t>  is  assumed.  Therefore, 
the  strength  degradation  scheme  will  look  like  the  following 

si=a+a-w,)5,  (3.io4) 

The  specific  details  of  property  degradation  scheme  used  in  this  work  including  the 
degradation  factors  used  for  the  different  materials  will  be  given  in  Section  9. 

The  two  degradation  schemes  involved  with  the  coupled  analysis,  that  is,  one  based  on 
the  stress  state/mechanical  damage,  and  the  other  based  on  the  oxidation,  need  to  be 
aggregated  to  provide  the  overall  mechanical  properties  of  the  material  based  on  the 
oxidation  level  and  the  mechanical  damage.  At  each  time  step,  this  overall  set  of 
properties  will  be  used  to  perform  the  stress  analysis  in  the  damage  progression  model, 
and  then  check  for  new  damage  based  on  the  failure  criteria.  In  this  work,  a  procedure 
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has  been  implemented  to  combine  the  two  degradation  schemes.  This  procedure  in  the 
coupled  analysis  algorithm  would  correspond  to  the  box  in  Figure  3.6  that  is  labeled 
‘Modify  mechanical  properties  based  on  current  oxidation  state  and  damage  state’.  Let 
us  consider  the  procedure  for  a  material  point  in  the  configuration.  The  procedure  takes 
the  initial  mechanical  properties  for  the  material  and  the  current  oxidation  and  damage 
state  as  input  and  returns  the  modified  mechanical  properties.  The  procedure  is  as 
follows.  The  initial  mechanical  properties  are  modified  following  the  degradation 
scheme  based  on  the  oxidation  level.  At  the  end  of  this  first  step,  the  compliance  matrix 
has  been  modified  according  to  Eq..(3.102),  and  the  strengths  have  been  modified 
according  to  Eq.(3.104)  In  the  second  step  the  new  properties  are  then  modified  again 
based  on  the  degradation  scheme  based  on  mechanical  damage.  Therefore,  at  the  end  of 
the  second  and  final  step,  the  properties  obtained  from  the  first  step  are  then  modified 
according  to  Eq.(3.6).  In  reality,  the  order  of  the  steps  do  not  matter  and  the  overall 
elastic  moduli  can  be  summarized  as  follows 


En  = 

"(1  +  0-0)^)" 
ax 

En> 

E22  — 

~(l  +  (l-0)/72)' 

a2 

E '22  ’  ^ 33 

"(l  +  a-O)^)" 

a3 

G12  = 

'(i+(i-o  )bA) 
a4 

g12. 

g23  = 
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a5 
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The  overall  strengths  would  be  simply  those  given  by  Eq..(3.104)  because  the 
degradation  scheme  based  on  mechanical  damage  does  not  modify  the  strengths  of  the 
material. 


41 


Figure  3.6:  Algorithm  for  one-way  coupled  oxidation-damage  progression  analysis 
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3.7.  Summary 

The  common  aspects  of  some  of  diffusion  and  oxidation  analyses  were  discussed  in  this 
section  and  a  common  analysis  procedure  was  determined  that  can  be  used  to  help 
design  an  analysis  framework.  The  procedure  is  then  used  to  derive  the  theory  and 
equations  involved  in  the  different  analysis  models  used  in  this  work.  The  finite  element 
formulations  for  the  models  were  derived  and  the  algorithms  for  the  analysis  were 
discussed.  It  also  discusses  the  strategies  involved  in  analyzing  periodic  configurations 
and  obtaining  effective  properties  for  periodic  microstructures.  The  models  described  in 
this  section  are  implemented  in  a  finite  element  analysis  framework  as  described  in  the 
next  section. 
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4.  Implementation 

This  section  describes  the  implementation  of  the  finite  element  formulation  presented  in 
Section  3.  The  overall  challenges  of  implementing  an  oxidation  analysis  into  a  finite 
element  code  will  be  discussed.  An  overall  analysis  algorithm  will  be  presented, 
followed  by  more  detailed  algorithms  for  various  subroutines  including  assembly, 
calculation  of  the  stiffness  matrix  and  load  vector,  and  post-processing  to  update 
oxidation  state. 

4.1.  Challenges 

The  finite  element  modeling  of  oxidation  is  a  transient  procedure  that  requires  the 
modeling  the  mechanisms  of  diffusion  of  oxygen  through  a  material  as  well  as  oxidation 
which  acts  as  a  sink  in  which  reaction  consumes  oxygen.  Modeling  diffusivity  is 
relatively  straightforward  with  the  main  challenges  pertaining  to  oxidation.  These 
challenges  include: 

•  Modeling  the  oxidation  process  (reaction  rate)  at  a  point 

•  Modifying  the  diffusivity  due  to  changes  caused  by  oxidation 

•  Developing  a  metric  (oxidation  state)  to  characterize  how  much  oxidation  has 
occurred  at  a  point 


4.2.  Overall  analysis  algorithm 

This  subsection  gives  the  basic  algorithm  that  drives  the  oxidation  analysis.  For 
simplicity  this  algorithm  is  for  manual  time  stepping  using  predefined  time  step  sizes. 
There  is  no  check  for  the  suitability  of  a  particular  time  step  to  give  an  acceptable 
change  in  the  solution.  For  efficiency  concerns,  it  is  likely  that  an  adaptive  time  step 
procedure  will  need  to  be  used. 

For  a  particular  time  step,  the  time  integration  parameters  c/,  and  a2  must  be  calculated 
using  Eqs.  (3.33)  and  (3.34)).  Next  the  equations  are  set  up  by  assembling  the  elemental 
M  and  F  (calculated  using  Eqs.  (3.97)  and  (3.98))  into  the  global  system  of  equations  to 
which  boundary  conditions  are  applied.  The  system  of  equations  is  solved,  with  the 
solution  being  the  increment  in  the  solution  between  the  current  time  step  and  the 
previous  time  step.  From  this  change  in  solution,  the  current  solution  is  calculated.  Using 
the  current  solution,  the  oxidation  state  is  updated  which  will  be  used  in  post-processing 
as  well  as  the  calculation  of  diffusivities  and  reaction  rates. 
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FOR  s  =  1  to  #  of  time  steps 

-  Calculate  time  integration  parameters 

-  Set  up  equations 

Assemble  M  and  F 
Apply  boundary  conditions 

-  Solve  M A q  =  F 

-  qs+ i  =  qs+hi 

-  Update  oxidation  state 

-  Output  (if  desired) 

END  FOR 


4.3.  Assembly 

This  subsection  gives  the  algorithm  that  assembles  the  equations  for  the  finite  element 
implementation  of  the  oxidation  analysis.  For  a  particular  element,  the  degree  of 
freedom  (DOF)  list  which  relates  the  local  DOF  numbering  of  the  element  to  the 
corresponding  global  DOF  is  calculated.  The  element  K  matrix  is  calculated  using  Eq. 
(3.77)  .  This  calculation  is  dependent  on  diffusivity  which  is  now  updated  to  reflect  the 
most  recently  available  oxidation  state  at  the  time  the  K  matrix  is  calculated.  The 
element  M  matrix  is  calculated  using  Eq.  (3.76)  which  is  independent  of  material 
properties  and  remains  constant  (for  each  element)  throughout  the  analysis.  The  element 
reaction  vector  R*  is  calculated  using  Eq.  (3.78)  and  is  dependent  on  temperature  and 
oxygen  concentration.  The  element  K,  M,  and  R*  are  utilized  to  form  M  and  F  as 
shown  in  Eqs.  (3.97)  and  (3.98).  The  degree  of  freedom  list  is  then  used  to  assemble  the 
element  matrix  and  vector  in  the  global  system  of  equations. 

FOR  i  =  1  to  #  of  elements 

-  Calculate  DOF  List 

-  Calculate  Element  K 

-  Calculate  Element  M 

-  Calculate  Element  R* 

-  Form  element  M  and  F 

-  Assemble  M  and  F 

END  FOR 
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4.4.  Stiffness  matrix  calculation 

This  subsection  details  the  calculation  of  the  elemental  K  matrix  (Eq.  (3.77)).  As  typical 
with  finite  element  formulations,  shape  function  derivatives  and  volumes  associated  with 
the  an  integration  point  must  be  calculated.  In  particular,  the  diffusivity  of  a  particular 
integration  point  must  be  calculated.  This  is  done  by  scaling  between  unoxidized  and 
oxidized  diffusivities  using  the  oxidation  state  at  the  integration  point  (see  Eq.  (4.1)). 
Using  this  and  the  aforementioned  calculation,  the  particular  integrations  points 
contribution  to  the  elemental  K  matrix  may  be  calculated. 


Dv  = 


1  -4 


D  + 
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(4.1) 


FOR  j  =  1  to  #  of  integration  points 

-  Calculate  volume  associated  with  integration  point 

-  Calculate  shape  function  derivatives 

-  Calculate  integration  point  diffusivity 

-  Calculate  integration  point's  contribution  to  element  K 

END  FOR 


4.5.  Reaction  vector  calculation 

This  subsection  describes  the  elemental  calculation  of  the  reaction  vector  R*  (Eq.  (3.78) 
).  Again  the  volume  associated  with  the  integration  point  and  the  shape  function  values 
are  calculated.  The  most  recent  nodal  concentrations  are  then  interpolated  for  the 
concentration  at  the  integration  point.  This  concentration  is  used  to  evaluate  the  function 
/(C)  using  the  material  parameter  (5  (Eq.(3.56)).  This  function  is  then  used  with  the 

temperature  dependent  reaction  rate  R0 ,  and  the  oxidation  state  at  the  integration  point 

to  calculate  a  reaction  rate  (Eq.  (3.55))  which  reflects  the  most  recent  oxygen 
concentration,  oxidation  state,  and  temperature. 

FOR  j  =  1  to  #  of  integration  points 

-  Calculate  volume  associated  with  integration  point 

-  Calculate  shape  function  values 

-  Interpolate  element's  nodal  concentration  values  at  integration  point 

-  Evaluate /(C)  function 

-  Calculate  reaction  rate  with  oxidation  state, /(C),  and  Ra 
END  FOR 
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4.6.  Updating  oxidation  state 

This  section  describes  the  post-processing  step  of  updating  the  oxidation  state.  After  the 
solving  step  is  completed,  the  oxidation  state  must  be  updated  for  every  integration 
point  in  the  model.  This  requires  an  integration  point  loop  nested  inside  a  element  loop 
as  shown  in  the  algorithm  below.  For  a  particular  integration  point,  an  updated  reaction 
rate  will  be  calculated  which  incorporates  the  most  current  oxygen  concentration  and 
oxidation  state  (as  well  as  temperature)  for  the  integration  point  of  interest.  The 
proportionality  parameter  a  is  calculated  for  the  particular  time  as  described  in  Section 
3.4.8.  Next  the  oxidation  state  integral  is  updated  to  reflect  the  current  time  step  (Eq. 
(4.2))  and  the  oxidation  state  is  calculated  using  Eq.(3.62) . 

ts+\ 

J  a(OR(Od£  *a(t,+1)R(t,+l)Al 

ts 

FOR  i  =  1  to  #  of  elements 

FOR  j  =  1  to  #  of  integration  points 

-  Evaluate  integration  point's  reaction  rate 

-  Calculate  proportionality  parameter  at  current  time 

-  Evaluate  current  time  step's  contribution  to  oxidation  state 

integral 

-  Calculate  oxidation  state  (Eq.  3.62) 

END  FOR 
END  FOR 

4.7.  Summary 

Section  3  provided  a  detailed  finite  element  formulation  for  oxidation  analysis.  This 
section  has  provided  the  details  for  implementing  such  a  formulation  into  a  finite 
element  code.  The  challenges  of  such  an  implementation  were  addressed. 
Implementation  algorithms  were  provided  that  outlined  the  more  challenging  details  of 
the  implementation  that  were  unique  from  a  typical  finite  element  implementation.  With 
the  implementation  of  an  oxidation  finite  element  formulation  complete,  the  next  section 
will  focus  on  ways  to  expedite  an  oxidation  finite  element  analysis. 


(4.2) 
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5.  Strategies  to  Expedite  Analysis 


Restriction  on  element  and  time  step  size  can  make  diffusion  and  oxidation  analysis  time 
consuming.  Oxidation  analysis  requires  a  longer  analysis  time  due  to  the  reaction  of 
oxygen  with  the  polymer  which  serves  to  slow  the  progression  of  oxygen.  In  order  to 
expedite  the  analysis  a  variety  of  strategies  were  considered.  These  include  optimization 
of  element  size  and  time  step  size,  an  adaptive  meshing  strategy,  and  a  decoupled 
subdomain  strategy  for  textiles.  The  use  of  approximate  boundary  conditions  was  also 
considered  to  reduce  the  size  of  the  analysis  region.  An  adaptive  time  stepping  algorithm 
has  also  been  developed.  Finally,  the  efficiency  of  the  implementation  has  been 
increased  mainly  through  parallel  algorithms  and  the  exploitation  of  direct  solver 
characteristics. 

5.1.  Optimization  of  element  size  and  time  step  size 

As  with  typical  mechanical  analyses,  the  accuracy  of  the  solution  depends  on  several 
factors,  one  of  which  is  the  element  size.  In  the  case  of  transient  analyses  like  diffusion 
or  oxidation,  the  solution  also  depends  on  the  time  step  size.  Depending  on  the  material 
properties  and  other  values  in  the  finite  element  formulas,  there  are  limits  to  the  element 
size  and  time  step  size  beyond  which  meaningless  results  are  obtained.  Analyses  were 
performed  using  COMSOL  Multiphysics  to  confirm  that  other  finite  element  packages 
had  the  same  limitations.  In  addition  to  the  basic  approximation  for  the  time  integrations, 
there  are  several  approximations  made  in  the  finite  element  formulation  to  handle  the 
nonlinearity  in  the  governing  equations.  The  accuracy  of  these  approximations  depends 
on  parameters  such  as  the  time  step  size  as  well. 

In  general,  the  optimal  time  step  size  need  not  be  constant  throughout  the  simulation 
because  of  the  nonlinear  oxidation  behavior.  This  means  that  the  time  step  size  can 
potentially  be  ramped  up  or  down  as  the  simulation  is  in  progress  so  as  to  maintain  the 
optimal  time  step  size.  To  summarize,  the  following  optimizations  can  be  made  to  an 
oxidation  simulation  in  order  to  make  it  run  more  efficiently: 

1.  Optimal  element  size 

2.  Optimal  time  step  size 

3.  Optimal  time  step  size  ramping 

Parametric  studies  were  conducted  to  determine  the  optimized  parameters  for  the 
materials  that  would  be  analyzed  in  this  work.  The  latter  part  of  this  section  will  discuss 
the  results  of  these  parametric  studies. 

The  remainder  of  this  section  describes  the  oxidation  behavior  in  neat  PMR-15  resin. 
Certain  characteristics  of  the  oxidation  behavior  can  be  exploited  to  develop  a  strategy  to 
speed  up  the  analysis.  For  this  purpose,  oxidation  of  a  simple  configuration  is 
considered.  The  simple  configuration  is  a  block  of  neat  resin  that  is  exposed  to  oxygen 
on  one  pair  of  opposite  surfaces  that  are  40  mm  apart  and  protected  from  oxygen  on  the 
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other  surfaces.  This  configuration  can  be  analyzed  using  a  1-D  model.  Moreover,  taking 
advantage  of  symmetry,  only  half  of  the  block  needs  to  be  modeled.  Table  5.1  gives  the 
material  properties  used  to  model  the  neat  PMR-15  resin.  For  a  complete  description  of 
the  different  oxidation  material  properties,  refer  to  Section  3.4.8.  A  uniform  element  size 
of  lmicron  and  time  step  size  of  0.15  minute  was  used  for  the  simulation. 

Figure  5.1  shows  the  predicted  oxidation  layer  growth  for  the  configuration  over  a 
period  of  200  hours.  Section  3.4.6  describes  how  the  oxidation  layer  growth  is 
determined.  It  can  be  seen  that  the  resin  oxidizes  very  quickly  in  the  initial  20  hours  or 
so  and  then  gradually  slows  down  to  where  the  oxidation  layer  grows  almost  linearly. 
Also  note  that  the  thickness  of  zone  II  or  the  active  zone  remains  fairly  constant 
throughout  the  entire  process. 


Table  5.1:  Oxidation  material  properties  for  neat  PMR-15  resin 


Neat  PMR-15  resin 

Diffusivity 

D unox 

Dox 

53.6xl0"6  mm2/min 

7 8. 22x1 O'6  mm2/min 

Ro 

3.5  mol/(m3min) 

<fioX 

0.187 

c 

0.79  mol/m3 

cc 

0.01-0.0067(t/40) :  t  <  40 
0.0033  :  t  >  40  (t  in  hours) 

m 

2 pC  [  pc  1 

1  +  /?c[  2{\  +  pC)\ 

P 

0.919 

The  difference  between  oxidation  and  diffusion-only  is  that  for  oxidation,  the  oxygen 
molecules  do  not  diffuse  as  quickly  because  they  are  consumed  in  oxidizing  the  material. 
Thus,  the  reaction  term  in  the  governing  equations  gives  the  effect  of  a  ‘moving  barrier’ 
that  allows  almost  no  oxygen  to  cross  over  to  the  other  side  of  the  active  zone  until  there 
is  a  sufficient  level  of  oxidation  within  the  active  zone.  This  is  evident  by  looking  at  the 
concentration  profiles  across  the  model  at  different  snap  shots  during  the  simulation. 
Figure  5.2  shows  the  concentration  profiles  in  the  model  at  t=2.5  hrs,  50  hrs  and  100  hrs. 
It  can  be  seen  that  all  the  profiles  have  a  similar  shape.  The  profiles  drop  almost  linearly 
from  the  exposed  edge  up  to  the  ‘moving  barrier’  and  the  concentration  is  practically 
zero  for  the  rest  of  the  model.  The  difference  in  each  profile  is  that  as  time  passes,  the 
location  of  the  ‘moving  barrier’  shifts  in  the  direction  of  the  oxygen  flow.  This 
movement  of  the  barrier  is  very  slow  compared  to  the  diffusion-only  process.  This  is 
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illustrated  in  Figure  5.2  by  the  concentration  profile  of  the  corresponding  diffusion 
model  at  15  minutes.  It  shows  that  with  only  15  minutes  of  diffusion,  the  oxygen 
concentration  at  every  point  in  the  model  has  already  surpassed  that  of  the  oxidation 
model  at  2.5  hours.  Even  after  100  hours  of  oxidation,  the  oxygen  concentration  is  still 
practically  zero  past  0.06  mm  whereas  the  corresponding  concentration  from  the 
diffusion  model  after  15  minutes  is  more  than  0.025  at  0.06  mm.  This  also  explains  why 
there  is  a  close  to  linear  drop  of  the  concentration  from  the  exposed  edge  to  the  ‘moving 
barrier’.  In  each  snapshot  of  concentration  profile  in  the  oxidation  process,  the  region  to 
the  left  of  the  moving  barrier  can  be  considered  as  a  diffusion  only  region  with  fixed 
concentration  boundary  conditions  -  the  specified  concentration  at  the  exposed  boundary 
and  zero  concentration  at  the  location  of  the  barrier.  Since  the  barrier  is  moving  very 
slowly,  the  concentration  profiles  at  the  various  time  steps  look  very  similar  to  that  for 
the  corresponding  diffusion-only  problem  at  steady-state,  which  is  a  nearly  linear 
variation  of  the  concentration.  Examination  of  this  behavior  gave  way  to  a  strategy  to 
further  expedite  the  oxidation  simulation.  This  strategy  was  called  the  "Adaptive 
Meshing  Strategy"  and  is  described  in  detail  in  the  next  section. 


Figure  5.1:  Predicted  oxidation  layer  growth  (Zone  I+II,  Zone  II)  in  neat  PMR-15 

resin 
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Figure  5.2:  Concentration  profiles  for  oxidation  and  diffusion  models. 


5.2.  Adaptive  meshing  strategy 

The  fact  that  the  concentration  of  oxygen  in  the  un-oxidized  region  of  the  material  is 
practically  zero  can  be  exploited  to  speed  up  the  analysis  by  constraining  the  degrees  of 
freedom(DOF)  in  most  of  the  un-oxidized  region  to  zero.  This  can  lead  to  a  considerable 
reduction  in  the  number  of  unknowns  to  be  solved  for,  especially  in  the  initial  period  of 
oxidation  because  most  of  the  material  is  un-oxidized  at  that  time.  The  challenge  is  in 
determining  which  regions  of  the  material  should  be  constrained  and  developing  an 
efficient  algorithm  so  that  this  can  be  automated.  The  regions  very  close  to  the  active 
zone  should  not  be  constrained  since  the  active  zone  is  slowly  moving  to  the  interior  of 
the  material  with  each  time  step  and  that  can  affect  the  solution.  Also,  the  regions  should 
not  be  permanently  constrained  because  that  implies  that  those  regions  will  never  get 
oxidized,  which  is  not  the  case. 

Based  on  these  requirements,  the  following  algorithm  was  developed  to  automatically 
determine  the  regions  to  be  constrained.  A  very  small  concentration  value  close  to  zero 
is  chosen,  say  C°,  in  order  to  determine  which  regions  are  to  be  constrained.  If  the 
concentration  at  a  node  is  more  than  C° ,  then  that  location  is  assumed  to  be  inside  the 
oxidation  layer  or  close  to  it  and  therefore  the  DOF  for  that  node  is  left  unconstrained. 
On  the  other  hand,  if  the  concentration  at  a  node  is  less  than  C°,  then  the  node  is 
assumed  to  be  in  the  un-oxidized  region  and  far  enough  from  the  active  zone,  therefore 
that  DOF  is  constrained.  This  check  is  not  performed  at  every  time  step.  Instead,  the 
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check  is  performed  every  15  or  20  time  steps  or  some  optimal  number  of  time  steps  (say, 
N)  chosen  depending  on  the  rate  the  active  zone  is  moving.  Therefore,  once  a  check  is 
performed,  the  constrained  DOFs  remain  constrained  for  the  subsequent  time  steps  until 
the  time  step  right  before  the  next  check.  In  this  time  step  preceding  the  check,  all  the 
artificial  constraints  are  removed  and  the  full  system  of  equations  is  solved.  This  allows 
a  minute  amount  of  oxygen  to  enter  the  previously  constrained  region.  In  the  next  time 
step,  the  check  is  performed,  at  which  time  some  of  the  previously  constrained  DOFs 
will  be  unconstrained  because  the  oxygen  concentration  has  increased  by  a  small 
amount.  This  cycle  is  repeated  throughout  the  simulation.  This  strategy  speeds  up  the 
analysis  by  a  large  factor  because  in  the  standard  analysis,  every  time  step  involves 
solution  of  the  entire  system  of  equations  whereas  in  the  adaptive  mesh  analysis,  the 
entire  system  of  equations  is  solved  only  every  N  time  steps.  During  the  other  time  steps, 
the  system  of  equations  solved  is  much  smaller.  The  check  to  determine  the  region  to  be 
constrained  is  also  performed  only  every  N  time  steps  and  the  computation  effort  used 
for  the  check  is  miniscule  compared  to  the  savings  obtained  by  solving  a  smaller  set  of 
equations.  In  addition  to  those  savings,  whenever  the  check  is  performed  and  a  region  of 
the  un-oxidized  material  is  constrained,  the  corresponding  elements  are  also  deactivated 
thereby  speeding  up  the  finite  element  assembly  process  as  well. 

The  choice  of  the  value  of  C°  has  an  effect  on  the  analysis  because  if  the  value  is  too 
large,  regions  that  are  close  to  the  active  zone  will  be  constrained  whereas  if  the  value  is 
too  small,  a  smaller  region  is  constrained  and  the  strategy  is  not  used  to  its  maximum 
potential.  Similarly,  the  number  of  time  steps  that  is  skipped  before  a  check,  N,  also  has 
an  effect  on  the  efficiency  of  the  simulation.  Parametric  studies  were  performed  by 
varying  the  two  parameters,  C°and  N  on  1-,  2-  and  3-D  models.  The  results  of  this 
parametric  study  are  presented  in  the  validation  section  of  this  report. 


5.3.  Decoupled  subdomain  strategy 

Textile  composites  have  multiple  microstructural  scales  -  the  fiber/matrix  scale,  the  tow 
architecture  scale  and  laminate  scale.  As  mentioned  in  the  previous  sections,  it  is  not 
practical  to  discretely  model  all  the  fibers  in  the  composite  because  of  modeling  and 
computational  challenges.  Effective  oxidation  material  properties  that  are  derived  in 
Section  3  and  validated  in  Section  7  are  used  to  model  the  tows  in  the  textile  composite. 
The  adaptive  meshing  strategy  described  in  Section  5.2  gives  considerable  savings 
compared  to  the  standard  finite  element  method  but  unfortunately,  it  is  not  enough  to 
make  the  3-D  analysis  of  textile  composites  feasible.  Given  the  length  scales  involved 
and  the  limitations  on  the  element  size,  the  mesh  would  require  a  huge  number  of 
elements.  This  would  make  generating  the  models  extremely  challenging,  and  analyzing 
the  models  practically  impossible.  Moreover,  considering  that  the  overall  goal  of  this 
research  effort  is  to  couple  the  oxidation  analysis  with  the  damage  progression  analysis, 
the  combination  would  be  prohibitively  expensive.  In  an  effort  to  make  this  more 
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feasible,  a  decoupled  subdomain  strategy  was  developed  to  make  the  oxidation  analysis 
more  efficient. 

The  strategy  applies  to  composite  laminates  that  are  exposed  to  oxygen  from  the  top  or 
bottom  (or  both)  surfaces,  but  not  the  lateral  surfaces.  The  strategy  is  illustrated  in 
Figure  5.3,  which  shows  a  l/8th  unit  cell  of  a  plain  weave  symmetrically  stacked  2-ply 
laminate.  The  decoupled  subdomain  strategy  takes  the  three-dimensional  model  and 
divides  it  up  into  individual  analysis  domains  in  the  in-plane  dimensions  as  shown  in 
Figure  5.3.  The  strategy  assumes  that  because  of  the  boundary  conditions  applied  on  the 
model,  the  oxidation  behavior  will  be  such  that  the  neighboring  domains  do  not  have  an 
effect  on  each  other,  essentially  assuming  that  oxygen  does  not  flow  from  one  domain  to 
another.  Therefore,  the  individual  domains  can  be  analyzed  separately.  Each  individual 
domain  is  a  three  dimensional  heterogeneous  analysis  region  with  curved  material 
boundaries  because  of  the  undulation  of  the  tows  in  the  textile  composite. 

The  strategy  assumes  that  the  change  in  the  diffusivity  due  to  the  undulation  is  not 
significant  because  the  rotation  angles  in  actual  composites  are  relatively  small.  The 
analysis  also  assumes  that  the  undulations  of  the  tows  are  not  significant  enough  to  cause 
an  impact  on  the  oxidation  behavior.  This  assumption  has  been  validated  and  is 
discussed  in  the  next  section.  Based  on  this  assumption,  the  individual  3-D  domain  can 
be  converted  into  an  equivalent  domain  with  straight  horizontal  material  boundaries 
based  on  the  volume  fraction  of  the  different  constituents  in  the  domain  as  illustrated  in 
Figure  5.3.  Since  the  new  equivalent  domain  has  no  inclined  material  boundaries,  it  can 
be  analyzed  with  a  simple  1-D  model.  Thus,  the  3-D  model  shown  in  Figure  5.3  can  be 
replaced  by  an  array  of  64  1-D  models,  thereby  reducing  analysis  time  significantly.  The 
decoupled  subdomain  strategy  is  implemented  in  the  finite  element  analysis  package  in 
such  a  way  that  the  input  to  the  model  is  the  same  as  the  conventional  3-D  model. 
Additional  pre-processing  work  is  not  required  and  the  array  of  1-D  models  is 
automatically  generated  and  analyzed  without  the  need  for  human  interaction.  Moreover, 
the  1-D  models  can  be  run  in  parallel  on  multi-core  processors,  thereby  increasing  the 
efficiency  even  further.  This  modeling  strategy  was  validated  by  using  a  2-D 
configuration.  The  validation  including  discussion  of  some  of  the  oxidation  behavior  is 
described  in  the  next  section. 


53 


3-D  analysis 


Equivalent  1-D 
domain 


Figure  5.3:  Schematic  of  hybrid  model  for  analyzing  textile  composites 

5.4.  Use  of  approximate  boundary  conditions 

Restrictions  on  element  size,  time  step  size,  and  longer  analysis  times  proved  to  be 
computationally  demanding.  The  unit  cell  for  a  symmetrically  stacked  plain  weave  is 
shown  in  Figure  5.4.  In  the  case  that  the  top  and  bottom  surfaces  are  exposed  to  an 
oxidizing  environment,  one  may  exploit  symmetries  and  reduce  the  size  of  the  analysis 
region  to  1/8  of  a  unit  cell  as  shown  in  Figure  5.5.  This  reduction  does  not  violate  any  of 
the  boundary  conditions  imposed  on  the  unit  cell  since  symmetry  is  being  exploited. 
Thus  it  is  said  to  have  "exact  boundary  conditions".  However,  element  size  restrictions 
combined  with  complex  weave  architecture  still  resulted  in  a  highly  refined  mesh  that 
made  timely  analysis  prohibitive.  Therefore,  the  analysis  region  was  further  reduced  to 
1/32  of  the  original  unit  cell  size  (see  Figure  5.6).  This  further  reduction  significantly 
decreases  the  analysis  region  but  does  not  allow  for  cross-flow  of  oxygen  between 
regions  as  the  1/8  model  does.  Thus  it  is  said  to  have  "approximate  boundary 
conditions".  However,  this  may  be  an  acceptable  approximation  if  the  cross-flow  for 
particular  configurations  is  minimal.  Validation  tests  will  assess  the  severity  of  this 
approximation.  Figure  5.7  shows  the  typical  mesh  refinement  that  is  used  in  an  oxidation 
analysis  which  results  in  a  relatively  high  number  of  degrees  of  freedom  given  the 
transient  nature  of  the  analysis. 
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Figure  5.7:  Typical  mesh  refinement  for 
1/32  unit  cell  (36000  DOF) 


Figure  5.6:  1/32  unit  cell  compatible  with 

approximate  boundary  conditions 


5.5.  Adaptive  time  stepping 

Manual  time  step  specification  and  time  step  ramping  allows  for  particular  segments  of  a 
transient  analysis  to  be  performed  with  varying  time  step  sizes  which  may  be  adjusted  to 
maintain  a  desired  level  of  accuracy  (compared  to  a  reference  solution  with  relatively 
small  time  steps)  while  increasing  the  overall  efficiency  of  an  analysis.  However,  the 
drawback  to  manual  time  stepping  is  that  previous  knowledge  of  the  transient  behavior 
of  a  configuration  is  typically  required.  Furthermore,  the  optimal  time  stepping 
determined  for  one  configuration  may  be  different  for  configurations  with  different 
geometry  or  constituents. 

Since  manual  time  stepping  may  essentially  require  parametric  studies  to  be  performed 
to  choose  optimal  parameters  for  a  given  configuration,  an  alternative  approach  is 
preferable.  An  adaptive  time  stepping  approach  that  actively  modifies  time  step  size  to 
maintain  a  maximum  value  of  change  in  solution  per  time  step  can  increase  efficiency 
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while  maintaining  accuracy.  The  proposed  adaptive  time  stepping  algorithm  ensures  that 
the  maximum  change  in  any  nodal  value  of  the  primary  variable  does  not  change  by 
more  than  a  user  specified  tolerance  during  a  single  time  step.  If  a  selected  time  step 
does  not  fulfill  the  requirement  imposed  by  the  tolerance,  the  time  step  is  rejected  and  a 
new  time  step  is  selected  based  off  of  Eq.(5.1).  This  procedure  is  repeated  until  a 
successful  time  step  is  executed  at  which  point  the  analysis  continues. 

The  new  time  step  is  computed  using  Eq.(5.1)  by  multiplying  the  old  time  step  by  the 
ratio  of  the  specified  TOL  and  the  maximum  absolute  change  in  the  primary  variable. 
The  factor  y  is  used  to  increase  the  efficiency  of  the  analysis  and  is  set  to  0.7  in  the 
current  implementation  to  slightly  reduce  the  new  time  step  with  the  goal  of  preventing 
unsuccessful  time  steps. 


At 


r  a*, 


TOL 


old 


max(l  Aq  I) 


(5.1) 


5.6.  Efficiency  improvements 

Even  with  optimization  of  element  size  and  time  step  size  through  ramping  or  adaptive 
time  stepping,  the  oxidation  analysis  of  all  but  the  simplest  of  configurations  proved  to 
be  quite  computationally  intensive.  Therefore,  further  increases  in  efficiency  were 
sought.  To  this  end,  the  code  was  re-examined  for  inefficiencies.  Many  modifications 
were  made,  but  the  most  significant  improvement  came  through  the  use  of  parallel 
algorithms  and  the  exploitation  of  some  characteristics  of  direct  solvers. 

5.6.1.  Parallel  algorithms 

The  relatively  recent  increase  of  parallel  computing  resources  makes  the  use  of  parallel 
algorithms  advantageous  for  increases  in  efficiency.  The  solvers  which  are  available  for 
oxidation  analysis  are  parallel.  Since  the  solving  step  is  typically  the  most  significant 
portion  of  analysis  such  a  capability  provides  an  obvious  advantage  over  solving  a 
system  of  equations  in  serial.  The  assembly  of  equations  is  typically  the  second  most 
expensive  portion  of  a  finite  element  analysis.  Parallel  algorithms  were  developed  using 
OpenMP,  which  also  allow  for  the  use  of  multi-point  constraints. 

The  use  of  multi-point  constraints  creates  an  increased  coupling  in  the  system  of 
equations  for  a  finite  element  analysis.  Typically  coupling  of  equations  is  more  or  less 
governed  by  the  connectivity  of  a  mesh.  However,  multi-point  constraints  can  couple 
two  degrees  of  freedom  which  have  no  geometric  connectivity  from  the  mesh. 
Therefore,  the  typical  use  of  domain  decomposition  to  split  the  mesh  between  processors 
is  not  as  useful  when  multi-point  constraints  are  involved.  Therefore,  a  new  parallel 
algorithm  was  developed  that  would  actually  split  the  system  of  equations  among  the 
processors  instead  of  the  elements.  Figure  5.4  visualizes  the  assignment  of  equations  to 
four  processors.  While  there  may  be  some  overlap  in  element  calculations  among  the 
processors,  this  method  ensures  multi-point  constraints  will  remain  intact  during  parallel 
assembly. 
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Figure  5.4:  Assignment  of  equations  to  processors 

With  this  methodology  two  parallel  assembly  algorithms  were  developed.  The  first 
algorithm  is  referred  to  as  the  "Elemental  Storage  Algorithm"  (ESA).  This  algorithm 
splits  the  elements  evenly  among  processors  for  calculation  and  stores  the  elemental  data 
in  memory.  A  second  element  list  for  every  processor  lists  the  elements  required  to 
assemble  the  equations  assigned  to  that  processor.  There  will  typically  be  some  overlap 
in  the  element  list  for  assembly  among  the  processors.  However,  there  are  no  concerns 
regarding  memory  access  since  no  processor  will  write  to  another  processor's  equations 
an  no  processor  will  calculate  the  same  element  as  another  processor.  While  it  is  possible 
that  multiple  processors  could  simultaneously  access  the  stored  elemental  data  from  the 
calculation  phase,  there  is  still  no  concern  since  the  processors  will  only  be  reading  the 
same  block  of  memory  and  not  attempting  to  write  to  it.  The  elemental  storage  algorithm 
has  the  advantages  of  scaling  very  well  with  the  number  of  processors  since  there  is  no 
overlap  in  the  element  calculations.  Furthermore,  no  synchronization  of  the  processors  is 
required.  However,  storage  of  elemental  data  can  be  impractical  for  larger  problems  due 
to  memory  restrictions.  This  algorithm  is  best  suited  for  small  to  moderate  sized 
problems  such  as  progressive  damage  or  transient  analysis.  The  elemental  storage 
algorithm  is  currently  implemented  for  oxidation  analysis. 

The  second  algorithm  is  referred  to  as  the  "On-the-Fly  Algorithm"  (OF A).  As  the  name 
suggests,  the  calculation  of  elemental  data  is  performed  on  the  fly  as  needed  for 
assembly  and  is  not  stored.  Therefore,  the  same  element  could  have  calculations 
performed  multiple  times  among  the  various  processors.  The  possibility  of  multiple 
processors  attempting  to  calculate  the  same  element  must  be  prevented  by  some 
synchronization  steps  which  can  reduce  the  efficiency  of  the  algorithm.  For  small  and 
moderate  sized  problems  OFA  does  not  scale  as  well  as  ESA  with  respect  to  the  number 
of  processors.  However,  as  the  problem  size  becomes  larger  there  is  significantly  less 
overlap  among  processors  and  the  speed  up  approaches  that  of  ESA.  Another  advantage 
of  OFA  is  that  there  are  no  memory  concerns  for  larger  models. 
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5.6.2.  Exploitation  of  direct  solvers 

Direct  equation  solvers  are  commonly  used  for  moderately  sized  system  of  equations 
and  can  offer  good  performance  compared  to  an  iterative  solver  for  the  same  system  of 
equations.  For  significantly  larger  systems  (greater  than  1  million  unknowns),  memory 
limitations  often  prevent  the  use  of  direct  solvers,  and  iterative  solvers  can  typically 
solve  larger  systems  in  less  time  than  direct  solvers.  However,  if  the  problem  size 
permits  the  use  of  a  direct  solver,  there  are  characteristics  that  may  be  exploited  to 
further  increase  the  efficiency. 

A  direct  solver  typically  has  three  phases: 

1.  Reordering  of  equations 

2.  Factorization/decomposition  of  equations 

3.  Back  substitution  to  solve  the  system  of  equations. 

The  first  and  last  phases  are  typically  cheaper  than  the  second  phase.  Throughout  the 
oxidation  analysis,  the  left  hand  side  matrix  will  be  updated  and  change  numerous  times. 
Therefore  the  factorization  and  back  substitution  phases  will  be  necessary  for  every 
solve  step.  However,  the  reordering  phase  depends  more  on  the  structure  of  the  matrix 
than  the  actual  values  of  the  matrix.  Since  the  general  structure  and  sparsity  of  the  matrix 
is  dictated  by  mesh  connectivity  and  boundary  conditions  this  structure  will  remain 
constant  throughout  the  typical  analysis.  Therefore,  the  reordering  may  be  performed 
once  and  stored.  This  can  be  a  significant  portion  of  the  analysis  for  larger  systems. 

5.7.  Summary 

This  section  has  detailed  some  of  the  attempts  and  methodologies  utilized  to  increase  the 
efficiency  of  oxidation  analysis.  This  includes  concerns  regarding  the  optimal  element 
size  and  time  step  size  along  with  the  use  of  time  step  ramping  to  define  different  time 
step  sizes  over  different  segments  of  the  analysis.  An  adaptive  meshing  strategy  was  also 
proposed  that  would  essentially  decrease  the  size  of  the  system  of  equations  depending 
on  how  far  oxygen  has  diffused  into  the  analysis  region.  In  the  case  of  textile  composite 
oxidation  analysis,  a  decoupled  subdomain  strategy  was  developed  which  converted  a 
three-dimensional  model  into  multiple  one-dimensional  domains.  Use  of  approximate 
boundary  conditions  to  reduce  the  problem  size  and  expedite  analysis  was  considered. 
To  provide  another  time  stepping  alternative,  an  adaptive  time  stepping  algorithm  was 
proposed  which  would  actively  change  the  time  step  size  while  maintaining  a  minimum 
change  in  solution  per  time  step.  Finally,  efficiency  improvements  including  the  use  of 
parallel  algorithms  and  the  exploitation  of  direct  solvers  was  discussed. 
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6.  Validation  of  Strategies  to  Expedite  Analysis 

Section  5  discussed  strategies  to  expedite  oxidation  analysis.  This  section  discusses  the 
validation  that  was  performed  on  such  strategies.  Also,  as  discussed  in  Section  3 
homogenized  properties  will  be  utilized  in  textile  oxidation  analysis  due  to  the  difficulty 
and  expense  of  modeling  the  microstructure  of  textiles. 

6.1.  Optimal  element  size  and  time  step  size 

Parametric  studies  were  conducted  to  determine  the  optimal  element  size  and  time  step 
size  as  well  as  the  time  step  size  ramping.  In  order  to  make  comparisons,  parametric 
studies  were  also  conducted  on  corresponding  diffusion  models.  Optimal  element  size 
and  time  step  size  were  determined  by  analyzing  the  same  configuration  described  in 
Section  5.1.  The  diffusivity  of  the  material  for  this  parametric  study  was  assumed  to  be 
53.6xl0'6  mm2/min,  which  is  the  diffusivity  of  the  un-oxidized  PMR-15  resin.  It  is 
important  to  note  that  this  parametric  study  is  not  extensive  and  does  not  look  all  the 
possible  parameters.  Therefore,  the  results  from  this  parametric  study,  in  essence,  are 
valid  only  for  material  properties  and  other  model  parameters  used  in  the  study.  In  order 
to  analyze  other  material  systems,  it  would  be  advisable  to  determine  the  optimal 
parameters  for  that  specific  system. 

One-dimensional  models  were  analyzed  using  various  element  sizes  and  time  step  sizes. 
The  reference  solution  was  assumed  to  be  that  obtained  from  using  linear  elements  with 
a  size  of  1  micron  and  a  time  step  size  of  0.15  minutes.  The  variation  of  average 
concentration  in  the  model  with  time  was  compared  for  the  different  models.  It  was 
observed  that  the  effect  of  the  element  size  and  time  step  size  on  the  results  were 
independent  of  each  other.  The  element  size  was  kept  constant  at  1  micron  and  models 
were  analyzed  with  varying  time  step  sizes  and  it  was  found  that  the  time  step  size  could 
be  raised  to  over  10  minutes  before  any  noticeable  difference  in  the  results  were 
observed.  When  the  time  step  size  was  kept  constant  at  0.15  minutes,  the  element  size 
could  be  increased  to  at  least  40  microns  without  any  noticeable  change  in  the  results.  A 
model  with  an  element  size  of  40  microns  and  a  time  step  size  of  10  minutes  also  yielded 
the  same  behavior  as  the  reference  model.  This  behavior  was  seen  for  both  linear  and 
quadratic  elements.  In  some  instances,  the  nodal  concentrations  drop  below  zero  but  they 
are  still  considered  numerical  zeros  and  these  negative  concentrations  do  not  have  any 
significant  effect  on  the  results.  The  same  results  were  obtained  when  a  parametric  study 
was  conducted  on  two  dimensional  models  with  eight-node  quadratic  elements.  A 
parametric  study  was  also  conducted  to  determine  the  effect  of  diffusivity  on  the 
allowable  time  step  size.  As  expected,  when  the  diffusivity  is  increased,  the  oxygen 
takes  less  time  to  saturate  the  material  and  the  optimum  time  step  size  required  in  order 
to  get  a  converged  solution  becomes  smaller.  It  was  also  found  that  increasing  the 
element  size  while  keeping  the  time  step  size  and  diffusivity  constant,  eventually  results 
in  negative  nodal  concentrations. 

Similar  to  what  was  done  for  the  diffusion  analysis,  the  optimal  mesh  size  and  time  step 
size  were  determined  by  analyzing  the  configuration  described  in  Section  5.1  using  the 
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material  properties  in  Table  5.1.  The  reference  solution  was  assumed  to  be  that  obtained 
from  using  linear  elements  with  a  size  of  1  micron  and  a  time  step  size  of  0.15  minutes. 
The  oxidation  layer  thicknesses  were  calculated  for  all  the  models  and  compared  to 
determine  the  accuracy.  The  oxidation  layer  consists  of  a  fully  oxidized  layer  (Zone  I) 
and  the  active  reaction  layer  (Zone  II).  Although  the  Zone  II  layer  is  defined  by  having 
an  oxidation  level  in  between  0  and  1,  for  all  the  oxidation  models  described  in  this 
paper,  a  tolerance  of  1%  is  allowed  on  those  limits.  Therefore,  an  element  is  assumed  to 
have  started  oxidizing  and  is  in  Zone  II  if  the  oxidation  level  at  each  of  the  material 
integration  points  falls  within  a  lower  limit  of  0.01  and  an  upper  limit  of  0.99.  If  the 
oxidation  state  is  above  0.99,  the  element  is  assumed  to  be  un-oxidized  and  if  it  is  below 
0.01  it  is  assumed  to  be  fully  oxidized.  A  post-processing  routine  was  written  that 
calculated  the  growth  of  the  oxidation  layer  along  a  line  in  a  model.  This  involved 
extrapolating  the  oxidation  state  values  from  the  integration  points  to  the  nodal  points, 
averaging  the  extrapolated  values  at  a  node  if  the  node  shared  elements  of  the  same 
material  and  solving  for  the  location  on  the  prescribed  line  where  the  oxidation  level 
value  met  the  specified  upper  and  lower  limits. 

Figure  6.1  shows  the  effect  of  the  size  of  linear  elements  on  the  oxidation  layer  growth 
with  a  constant  time  step  size  of  0.15  minutes.  It  shows  that  the  models  using  4-micron 
and  8-micron  size  elements  closely  agree  with  the  model  using  1 -micron  elements 
whereas  the  model  using  12-micron  elements  over  predicts  the  thickness.  The  model 
with  8-micron  elements  shows  a  distinct  oscillation  in  the  curve.  This  is  believed  to  be 
caused  due  to  errors  from  extrapolation  of  the  oxidation  state  values  from  the  integration 
points  to  the  nodal  points.  Nevertheless,  it  can  be  seen  that  upper  bound  of  the  curve  is 
very  close  to  the  results  of  the  1 -micron  size  model.  The  model  with  4-micron  elements 
shows  slight  oscillations  as  well  but  it  is  able  to  predict  the  thickness  growth  very  well. 

The  effect  of  the  time  step  size  was  also  investigated  by  keeping  the  element  size 
constant  and  varying  the  time  step  size.  Figure  6.2  shows  the  oxidation  layer  growth  for 
different  models  when  the  element  size  is  kept  constant  at  2  microns  and  the  time  step 
size  varies  from  0.15  mins  to  0.8  mins.  It  can  be  seen  that  the  time  step  size  can  be 
doubled  from  0.15  mins  to  0.3  mins  without  any  perceivable  effect  on  the  results.  When 
the  time  step  size  is  raised  to  0.5  mins,  some  difference  can  be  seen  in  the  initial  part  of 
the  simulation  while  the  latter  part  still  predicts  the  oxidation  growth  fairly  well. 
Increasing  the  time  step  size  to  0.8  mins  affects  the  results  considerably  especially 
during  the  initial  part  of  the  simulation.  This  kind  of  behavior  for  the  effect  of  time  step 
size  on  the  predicted  oxidation  growth  was  seen  for  both  linear  and  corresponding 
quadratic  elements.  The  trends  also  show  that  the  time  step  size  is  more  critical  to  the 
initial  part  of  the  simulation  where  the  oxidation  growth  is  nonlinear.  For  many  of  these 
models,  the  nodal  concentrations  calculated  would  be  numerical  zeroes  that  go  below 
zero.  When  the  program  encounters  such  values,  they  are  converted  to  zero  so  that  it 
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Figure  6.1:  Effect  of  element  size  on  oxidation  layer  growth  (Zone  I+II)  for 
neat  resin  (using  linear  elements  and  time  step  size  of  0.15  mins) 


Figure  6.2:  Effect  of  time  step  size  on  oxidation  layer  growth  for 
neat  resin  (using  2  micron  linear  elements) 
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does  not  use  negative  concentrations  in  the  calculation  of  the  reaction  terms  and  the 
oxidation  state,  which  would  physically  mean  a  reversal  of  the  oxidation  process. 

In  order  to  speed  up  the  analysis,  the  behavior  of  the  model  when  the  time  step  size  is 
gradually  increased  was  investigated.  As  seen  from  the  results  of  the  previous  parametric 
study,  a  time  step  size  of  no  more  than  0.3  minutes  was  required  to  accurately  model  the 
initial  part  of  the  simulation  where  layer  growth  is  highly  nonlinear.  The  layer  growth 
behavior  becomes  close  to  linear  once  the  model  has  undergone  oxidation  for  40  hours, 
which  is  when  the  time-dependent  material  property,  a  changes  from  decreasing  linearly 
with  respect  to  time  to  a  constant  value  of  0.0033.  Based  on  this,  a  parametric  study  was 
conducted  where  the  models  used  a  time  step  size  of  0.3  mins  for  the  initial  40  hours  of 
the  simulation  and  for  the  other  160  hours,  the  different  models  used  different  time  step 
sizes.  The  reference  model  used  a  time  step  size  of  0.15  mins  for  the  entire  200  hours. 
All  the  models  used  elements  with  a  size  of  2  microns.  Figure  6.3  shows  that  when  the 
time  step  size  is  ramped  up  from  0.3  minute  to  1  minute,  the  predicted  oxidation  growth 
curve  is  barely  distinguishable  from  that  of  the  reference  model.  The  results  are  fairly 
reasonable  even  when  the  time  step  size  is  ramped  up  to  5  mins.  As  shown  in  Figure  6.3, 


Figure  6.3:  Oxidation  layer  growth  (Zone  I+II,  Zone  II)  for  neat  resin  (using  2 
micron  linear  elements  and  time  step  size  of  0.30  mins  for  the  first  40  hours 
and  different  ramped  time  step  sizes  thereafter) 


the  differences  in  the  curves  are  considerable  when  the  time  step  size  is  ramped  to  10 
mins.  Figure  6.3  show  that  the  models  also  predict  the  Zone  II  thickness  fairly  well. 
Ramping  up  the  time  step  size  tremendously  reduces  the  computational  time  required  for 
the  analysis  compared  to  using  a  constant  time  step  size  of  0.3  minutes.  A  constant  time 
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step  size  of  0.3  mins  used  for  simulating  200  hours  of  oxidation  takes  up  40,000  time 
steps  whereas  using  a  model  that  uses  0.3  mins  for  the  first  40  hours  and  5  mins  for  the 
remaining  160  hours  takes  up  only  9920  time  steps.  This  makes  a  computational  savings 
of  over  75%. 


6.2.  Adaptive  meshing  strategy 


Parametric  studies  were  performed  to  determine  the  optimal  parameters  for  the  Adaptive 
Meshing  Strategy  as  well  as  potential  computational  savings.  The  one-dimensional 
configuration  in  section  5.1  is  analyzed  using  the  Adaptive  Meshing  Strategy  described 
in  Section  4.2.  The  two  parameters  that  were  varied  were  C°  and  N.  All  the  models  in 
this  particular  parametric  study  use  1  micron  size  elements  and  time  step  size  ramping 
where  the  first  40  hours  use  0.3  minute  time  steps  and  the  remaining  160  hours  use  1 
minute  time  steps.  The  oxidation  layer  growth  from  the  different  models  is  compared 
with  a  reference  model  that  uses  the  standard  oxidation  analysis.  Figure  6.4  shows  the 
oxidation  layer  growth  for  models  that  have  a  constant  C°  of  0.01  and  three  different  N 
values  of  50,  100  and  200.  It  shows  that  for  N  values  of  50  and  100,  the  oxidation  layer 
growth  predicted  is  very  close  to  that  of  the  reference  model.  Even  for  the  model  with  an 
N  value  of  200,  it  is  seen  that  there  is  close  agreement  till  about  40  hours  after  which  the 
time  step  size  is  ramped  up  to  1  minute.  This  indicates  that  the  value  of  the  threshold 
concentration,  C°  is  too  high  and  that  the  oxidation  front  is  creeping  up  to  the 
constrained  region  and  the  active  region  of  the  mesh  is  not  re-evaluated  quickly  enough. 
That  is  why  for  lower  N  values  such  as  50  or  100,  the  prediction  of  oxidation  layer 
growth  is  much  better.  This  means  that  if  the  threshold  concentration,  C°  is  lowered, 
that  would  make  the  constrained  region  smaller  thereby  taking  it  longer  for  the  oxidation 
front  to  reach  the  region.  Therefore,  lowering  C°  should  allow  increasing  N  while 
maintaining  the  accuracy.  This  was  validated  by  analyzing  a  similar  set  of  models  as 
earlier  except  with  a  C°  value  of  le-3.  Figure  6.5  shows  that  N  value  of  200  does  a  very 
good  job  of  predicting  the  oxidation  layer  growth  whereas  when  C°  had  a  value  of  le-2, 
an  N  value  of  200  did  a  very  poor  job  of  predicting  the  layer  growth.  Moreover,  even  an 
N  value  of  300  does  a  good  job  and  it  is  only  when  it  is  increased  to  400  that  the 
accuracy  of  the  prediction  even  starts  to  deteriorate.  This  same  trend  was  seen  when  the 
C°  was  lowered  even  further  to  10  4  and  10"5.  On  the  other  hand,  when  C°  is  lowered,  the 
constrained  region  is  reduced  and  the  maximum  potential  of  the  Adaptive  Meshing 
Strategy  is  not  achieved.  Figure  6.6  shows  computational  time  savings  achieved  when 
the  Adaptive  Meshing  Strategy  compared  to  the  corresponding  standard  analysis  with 
time  step  size  ramping.  The  computation  time  savings  is  defined  by  using 


computational  time  savings  = 


^  time  taken  by  Adaptive  Meshing  Strategy  ^ 
time  taken  by  standard  analysis 
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Due  to  constraints  with  the  computational  resources,  it  was  not  possible  to  obtain 
accurate  timings  of  the  analysis  but  it  still  gives  a  good  sense  for  the  trends  in  the 
savings  achieved  when  the  value  of  C°  is  lowered.  As  illustrated  in  the  figure,  as  the 
value  of  C°  is  lowered  from  le-2  to  le-5,  the  computation  time  savings  decreased  from 
-68%  to  -58%. 

In  order  to  see  how  this  analysis  strategy  fares  when  the  dimensionality  of  the  model  is 
increased,  the  same  configurations  were  analyzed  using  2D  and  3D  models.  The  2D 
mesh  had  dimensions  of  200  x  10  elements  using  8-noded  2D  elements  of  size  1  micron. 
The  2D  model  had  a  total  of  6421  DOFs.  The  3D  model  had  dimensions  of  5  x  5  x  200 
elements  using  20-noded  brick  elements  of  size  1  micron.  The  3D  model  had  a  total  of 
26496  DOFs.  Again  constraints  on  the  computational  resources  prevented  accurate 
timings  of  the  analyses  but  it  did  give  the  same  kind  of  trend  for  all  the  models  analyzed. 
Figure  6.7  shows  the  computational  time  savings  achieved  when  C°  was  kept  at  a 
constant  value  of  le-3  and  the  value  of  N  has  been  varied  for  the  corresponding  ID,  2D 
and  3D  models.  The  results  were  not  conclusive  enough  to  determine  any  strong  trends. 
In  general,  it  was  seen  that  the  percentage  savings  reduced  for  the  2D  model  compared 
to  the  ID  model.  On  the  other  hand,  the  3D  models  generally  gave  a  better  percentage 
savings  compared  to  the  ID  models.  It  is  estimated  that  this  trend  is  due  to  the  nature  of 
the  system  of  equations  related  to  ID,  2D  and  3D  models.  For  all  the  analyses  performed 
in  this  work  other  than  this  parametric  study,  C°  was  chosen  to  be  0.0001  mol/m3  and  N 
was  chosen  to  be  20  time  steps. 


Figure  6.6:  Computational  time  savings  for  parametric  study  of  ID  expedited 
analysis  models  with  various  C°  and  N  values) 


65 


Figure  6.7:  Computational  time  savings  for  parametric  study  of  ID,  2D  and  3D 
expedited  analysis  models  with  C°=  10'3  and  various  N  values) 


6.3.  Decoupled  subdomain  strategy 

The  decoupled  subdomain  strategy  was  examined  to  assess  the  ability  of  the  method  to 
characterize  a  two-dimensional  or  three-dimensional  oxidation  analysis  through  the  use 
of  multiple  one-dimensional  domains.  At  the  time  the  method  was  developed,  a  two- 
dimensional  validation  was  performed  since  the  initial  efficiency  of  the  implementation 
prohibited  oxidation  analysis  of  three-dimensional  weave  models.  However,  recent 
advancements  in  the  efficiency  have  allowed  for  this  method  to  be  validated  against 
three-dimensional  results.  It  was  found  the  decoupled  subdomain  strategy  predicts  the 
results  of  an  oxidation  analysis  very  well  compared  to  a  three-dimensional  solution. 

6.3.1.  2D  Validation 

A  simple  2-D  heterogeneous  configuration  with  two  materials  was  chosen  where  the 
material  boundary  is  straight  but  at  an  angle  to  the  horizontal  edge  as  shown  in  Figure 
6.8(a).  The  bottom  edge  is  assumed  to  be  exposed  to  oxygen  whereas  the  other  three 
edges  are  assumed  to  be  impermeable.  The  configuration  has  the  dimensions  200 
microns  x  100  microns.  The  material  in  the  lower  region  is  assumed  to  be  neat  PMR-15 
resin  and  the  other  material  is  assumed  to  be  a  homogenized  graphite/PMR-15  tow  with 
a  fiber  fraction  of  55.6%.  The  2D  plane  in  the  configuration  is  assumed  to  be  the  plane 
perpendicular  to  the  fiber  axis  in  the  tow  and  therefore  only  the  transverse  diffusivities 
of  the  tow  will  be  used  in  the  2D  analysis.  The  material  properties  of  the  tow  are 
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calculated  using  the  formulas  described  in  Section  3.4.7.  The  material  properties  of  the 
resin  are  given  in  Table  5.1  and  that  of  the  homogenized  tow  using  the  aforementioned 
formulas  are  given  in  Table  6.1.  The  region  is  first  divided  into  two  domains  and 


Equivalent  1-D 
Domains 


(b)  Regions  represented  by 
the  1-D  domain  1  and  2 


Domain  1  Domain  2 

(c)  Regions  represented 
when  domain  1  and  2 
are  further  subdivided 


Domain  1-1  Domain  1-2  Domain  2-1  Domain  2-2 


Figure  6.8:  2-D  configuration  for  validating  hybrid  model 

converted  into  equivalent  1-D  models  as  shown  in  Figure  6.8(b).  To  compare  the 
oxidation  layer  growth  predicted  by  the  1-D  models  with  the  behavior  in  the  actual  2-D 
model,  the  oxidation  layer  growths  along  different  vertical  lines  (numbered  in  Figure 
6.8(a))  in  the  2-D  model  are  compared.  Figure  6.9(a)  plots  the  oxidation  growth  given  by 
the  equivalent  1-D  domain  1  model  along  with  that  along  lines  1,  3  and  5.  It  shows  that 
the  1-D  result  agrees  very  closely  with  that  of  line  3  and  not  so  much  with  that  of  lines  1 
and  5,  which  are  on  the  extreme  edges  of  domain  1.  Similar  trends  are  seen  in  Figure 
6.9(b),  which  shows  corresponding  plots  for  domain  2.  The  domains  are  then  further 
subdivided  into  domains  1-1,  1-2,  2-1  and  2-2  as  shown  in  Figure  6.8(c).  The 
corresponding  oxidation  growth  plots  for  domains  1-1,  2-1,  1-2,  and  2-2  are  shown  in 
Figure  6.9(c),  (d),  (e),  and  (f)  respectively.  As  expected,  these  results  show  that  the 
equivalent  ID  domain  models  perform  better  at  simulating  the  oxidation  layer  growth 
when  the  domain  size  is  reduced. 
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Table  6.1:  Oxidation  material  properties  for  the  homogenized  tow  (Vf=55.6%) 


Homogenized  tow 
(Vf=55.6%) 

Transverse 

Diffusivity 

Dox 

33.07  xlO"6  mm2/min 
48.27  xlO'6  mm2/min 

R0 

1.554  mol/(m3min) 

K 

0.639 

c 

0.3507  mol/m3 

GC 

0.01-0.0067(t/40) :  t  <  40 
0.0033  :  t  >  40  (t  in 
hours) 

m 

hO 

n 

— i 

n 

_ i 

l  +  pc'l  2(1+^C)J 

p 

0.919 
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Oxidation  layer  thickness  (mm)  Oxidation  layer  thickness  (mm)  Oxidation  layer  thickness  (mm) 


(a)  Domain  1  (b)  Domain  2 


(c)  Domain  1-1 


(d)  Domain  2- 1 


(e)  Domain  1  -2  (f)  Domain  2-2 


Figure  6.9:  Comparison  of  the  oxidation  layer  growth  from  the  different  1-D 
models  with  the  growth  in  the  2D  configuration 
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One  interesting  behavior  that  was  noticed  during  the  validation  was  that  when  simulating 
oxidation  of  a  heterogeneous  model  with  neat  matrix  and  tow,  the  predicted  oxidation 
growth  seems  counter-intuitive  when  compared  to  that  of  a  model  with  neat  matrix 
alone.  Consider  the  equivalent  1-D  configuration  for  domain  1  shown  in  Figure  6.10, 
which  is  a  heterogeneous  model  with  neat  resin  and  homogenized  tow.  Figure  6.11 
compares  the  predicted  oxidation  layer  growth  for  the  configuration  in  Figure  6.10  with 
that  of  a  pure  resin  model.  One  would  intuitively  expect  that  since  the  model  with  the 
tow  is  assumed  to  have  inert  and  impermeable  fibers,  this  would  slow  down  the 
oxidation  layer  growth  compared  to  a  neat  resin  model  that  has  no  fibers.  But  Figure 
6.11  shows  that  the  model  with  the  resin  and  tow  has  a  faster  oxidation  layer  growth. 


- ►x 

02— ► 

Resin 

Tow 

0.06mm 

^  W 

0.04mm 

Figure  6.10:  Equivalent  ID  configuration  for  domain  1 


Figure  6.11:  Comparison  of  oxidation  layer  growth  in  the  domain  1  (resin/tow) 

model  and  neat  resin  model 


On  further  investigation,  it  was  seen  that  a  number  of  factors  influence  this  behavior. 
The  tow  in  the  model  acts  like  a  pseudo-barrier  allowing  the  resin  to  saturate  with 
oxygen  much  faster  than  the  tow.  Until  the  oxidation  front  reaches  the  vicinity  of  the 
material  boundary,  both  the  models  behave  in  the  same  manner  because  the  tow  has  no 


70 


effect  on  the  matrix  that  is  being  oxidized  ahead  of  it.  But  once  the  tow  begins  to  oxidize 
as  well,  the  interface  conditions  regulate  the  flow  of  oxygen  from  the  matrix  into  the  tow 
and  free  oxygen  starts  to  build  up  in  the  matrix.  This  is  evidenced  in  Figure  6.12  which 
shows  the  oxygen  concentration  profile  in  the  model  at  100  hours.  Figure  6.12  shows 
that  the  resin  region  in  the  resin/tow  model  (from  0  to  0.06  mm)  has  more  oxygen  than 
the  same  region  in  the  neat  resin  model.  The  oxygen  in  the  tow  region  (from  0.06  to  0.1 
mm)  is  also  more  than  that  in  the  same  region  for  the  neat  resin  model.  This  could  be  due 
to  a  combination  of  factors.  First,  note  that  at  100  hours,  the  oxidation  front  has  crossed 
the  material  boundary  but  is  not  too  far  from  it.  The  material  boundary  is  at  0.06mm  and 
the  oxidation  front  at  100  hours  can  be  considered  to  be  around  0.08mm,  beyond  which 
the  oxygen  concentration  is  practically  zero.  Secondly,  the  tow  has  less  amount  of  resin 
that  can  be  oxidized  and  therefore  the  maximum  reaction  rate  is  also  less  than  that  of  the 
neat 


x  (mm) 


Figure  6.12:  Comparison  of  concentration  profile  in  the  resin/tow  model  and 
neat  resin  model  at  100  hours 

resin.  That  also  means  that  the  region  consumes  less  oxygen  (for  oxidation)  than  the  neat 
resin.  Since  the  reaction  rate  in  the  tow  is  less  than  that  in  the  neat  resin  and  the 
oxidation  front  is  fairly  close  to  the  material  boundary,  the  tow  region  between  the 
material  boundary  and  the  oxidation  front  also  starts  accumulating  more  oxygen  than  the 
corresponding  region  in  the  neat  resin.  Figure  6.13  gives  the  amount  of  free  oxygen  in 
the  model  throughout  the  simulation.  It  shows  that  until  about  60  hours,  the  resin/tow 
model  and  the  neat  resin  model  have  the  same  amount  of  free  oxygen,  but  after  60  hours 
the 
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Figure  6.13:  Comparison  of  amount  of  free  oxygen  in  the  resin/tow  model  and 
neat  resin  model 


resin/tow  model  builds  up  more  oxygen  in  its  material.  This  is  not  to  be  confused  with 
the  amount  of  oxygen  consumed  in  oxidizing  the  polymer  in  the  resin  and  tow  regions. 
The  neat  resin  model  is  expected  to  consume  more  oxygen  than  the  resin/tow  model 
because  it  has  more  material  that  can  be  oxidized  and  this  is  shown  in  Figure  6.14.  Once 
the  oxygen  starts  to  build  up  in  the  matrix,  it  becomes  fully  oxidized  more  quickly  and 
all  the  incoming  oxygen  is  directed  into  the  tow  region,  which  is  then  used  up  to  oxidize 
the  polymer  in  the  tow.  Also  note  that  an  oxidation  level  of  0.99  at  a  material  point  in  the 
neat  resin  region  indicates  that  1%  of  the  resin  in  the  material  has  oxidized.  On  the  other 
hand,  saying  that  1%  of  the  resin  in  a  material  point  in  the  tow  region  corresponds  to  an 
oxidation  level  defined  by  Eq.(3.99),  which  for  this  model  is  0.99556.  Figure  6.15  shows 
the  oxidation  level  profile  in  the  model  at  100  hours.  The  inset  plot  in  Figure  6.15  shows 
a  close  up  of  the  oxidation  state  of  the  two  models  between  0.08mm  and  0.095mm.  It 
shows  that  the  oxidation  level  in  the  resin/tow  model  dips  below  0.99556  at  about 
0.085mm  (at  location  A)  whereas  in  the  neat  resin  model,  it  dips  below  0.99  at  about 
0.077mm  (at  location  B).  This  snapshot  of  the  simulation  at  100  hours  shows  what  the 
oxidation  profile  in  the  two  models  looks  like  when  the  oxidation  layer  thickness  in  the 
resin/tow  model  is  larger  than  that  in  the  neat  resin  model.  Overall,  this  oxidation 
behavior  in  the  resin/tow  model  is  due  to  a  combination  of  factors  such  as  the  effective 
oxidation  properties  of  the  tow  as  well  as  the  diffusion  behavior  in  heterogeneous 
models  and  the  relatively  slow  movement  of  the  oxidation 
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Figure  6.14:  Comparison  of  amount  of  oxygen  consumed  in  the  resin/tow  model 

and  neat  resin  model 


Figure  6.15:  Comparison  of  oxidation  level  (d>)  profile  in  the  resin/tow  model 
and  neat  resin  model  at  100  hours 


front  into  the  interior  of  the  material.  It  would  also  depend  on  the  volume  fraction  of  the 
constituent  materials  and  the  configuration  of  the  materials  in  the  heterogeneous  model. 
Therefore,  the  location  of  the  material  boundary  in  the  configuration  also  has  an  impact 
on  the  oxidation  behavior. 
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6.3.2. 


3D  Validation  of  DSS  and  use  of  approximate  boundary  conditions 


Increases  in  efficiency  made  oxidation  analysis  of  three-dimensional  weave  models 
feasible.  However,  analysis  still  came  with  significant  computational  expense. 
Therefore,  an  approximate  boundary  condition  strategy  was  utilized  to  reduce  the 
problem  size  to  approximately  1/4  of  that  using  exact  boundary  conditions.  Figure  5.4 
shows  a  unit-cell  for  the  configuration  of  interest  for  oxidation  analysis  of  a  plain  weave 
2-ply  symmetric  laminate.  Oxygen  is  entering  from  the  top  and  bottom  surfaces  and  all 
other  surfaces  have  zero  flux  due  to  periodicity.  Figure  5.5  shows  the  analysis  region 
(1/8  unit  cell)  with  exact  boundary  conditions  once  symmetries  in  the  configuration  are 
exploited.  Here  oxygen  is  flowing  from  the  top  surface  and  all  other  surfaces  have  zero 
flux  due  to  symmetry.  Figure  5.6  shows  the  1/32  unit  cell  using  approximate  boundary 
conditions  assuming  all  surfaces  except  for  the  top  have  zero  flux.  Note  that  oxidation 
meshes  used  were  much  more  refined  than  that  shown  in  Figures  5.6  and  5.7.  This  was 
done  to  more  clearly  illustrate  tow  architecture.  Figure  5.8  shows  a  typical  mesh 
refinement  used  for  oxidation  analysis  of  a  three-dimensional  weave  model.  This  mesh 
contains  a  maximum  element  thickness  of  5  microns,  resulting  in  approximately  36000 
degrees  of  freedom  (DOF). 

With  results  available  for  a  three-dimensional  weave  model,  the  approximate  boundary 
condition  strategy,  as  well  as  the  previously  developed  decoupled  subdomain  strategy 
can  now  be  evaluated.  For  the  oxidation  analysis  of  textile  composites,  the  oxidation 
layer  growth  (essentially  how  deep  oxidation  is  occurring  in  the  composite  with  time)  is 
determined  along  various  lines  through  the  thickness.  The  oxidation  envelope  is  defined 
as  the  minimum  and  maximum  oxidation  layer  thickness  at  a  given  time.  Figure  6.16 
shows  the  comparison  of  oxidation  envelopes  for  a  2-ply  symmetric  plain  weave 
laminate. 

The  relative  difference  between  the  minimum  and  maximum  oxidation  layer  thickness 
that  remains  approximately  constant  with  respect  to  time  up  to  160  hours.  At  50,  100, 
and  150  hours  the  maximum  layer  thickness  increases  28-29%  over  the  minimum  layer 
thickness.  At  200  hours  the  maximum  layer  thickness  increases  37%  over  the  minimum 
layer  thickness.  However,  note  that  after  160  hours  there  is  a  sharp  increase  in  the 
maximum  oxidation  layer  thickness.  This  is  believed  to  be  due  to  mesh  refinement  and 
weave  architecture  in  a  certain  region.  More  refined  meshes  will  need  to  be  analyzed  to 
investigate  oxidation  behavior  in  this  region. 

The  results  for  exact  boundary  conditions,  approximate  boundary  conditions,  and  the 
decoupled  subdomain  strategy  agree  remarkably  well.  While  there  are  noticeable 
differences  in  the  distributions,  it  is  emphasized  that  these  differences  are  typically  less 
than  the  maximum  element  thickness  used  in  a  model.  Therefore,  more  mesh 
convergence  studies  should  be  conducted  before  trying  to  explain  the  relatively  minimal 
differences  in  oxidation  layer  growth  distributions.  The  three-dimensional  results  also 
serve  to  assess  the  severity  of  the  assumptions  in  the  decoupled  subdomain  strategy. 
Considering  all  of  the  uncertainties  in  this  complex  process,  the  highly  efficient 
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decoupled  subdomain  strategy  yields  results  that  should  be  considered  quite  acceptable 
for  the  current  configuration.  Also,  if  more  accurate  results  are  required,  the  decoupled 
subdomain  strategy  provides  a  quick  way  to  identify  the  most  important  cases  to  study 
further  with  fully  3-D  analysis. 


Figure  6.16:  Comparison  of  oxidation  envelopes  for  a  plain  weave  laminate  up  to  200 

hours 

6.4.  Mesh  convergence  study  for  3D  analysis 

With  efficiency  increases  allowing  three-dimensional  analysis  of  textile  composites,  a 
mesh  convergence  study  was  in  order.  While  3D  analysis  was  still  time  consuming,  a 
few  mesh  refinements  were  considered  to  study  the  effect  of  refinement  on  the  predicted 
behavior.  Figure  6.17  shows  the  meshes  considered  for  a  convergence  study.  These 
meshes  are  characterized  by  the  maximum  dimensions  of  an  element  in  the  mesh.  For 
example,  a  mesh  characterized  by  "10  x  15  pm"  will  have  elements  with  in-plane 
dimensions  of  15  x  15  pm  and  a  maximum  element  thickness  of  10  pm.  The  metric  used 
to  measure  convergence  was  the  oxidation  envelope  that  was  used  in  the  comparison  of 
exact  and  approximate  boundary  conditions  in  the  previous  section.  Figure  6.18  shows  a 
comparison  of  oxidation  envelopes  for  the  meshes  shown  in  Figure  6.17.  All  the  meshes 
show  approximately  the  same  response  to  oxidizing  environment  with  regards  to  the 
oxidation  envelope.  There  is  a  larger  deviation  for  the  coarsest  mesh  (10  x  15  pm) 
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beyond  150  hours,  but  this  deviation  is  relatively  minimal.  Thus,  the  mesh  shown  in 
Figure  6.14a  will  likely  be  sufficient  for  oxidation  analysis  of  a  plain  weave  composite. 


c.)  5  x  15  pm  mesh  d.)  2.5  x  15  pm  mesh 

Figure  6.17:  Various  mesh  refinements  considered  for  three-dimensional  oxidation 

analysis  of  a  plain  woven  composite 
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Figure  6.18:  Comparison  of  oxidation  envelopes  for  a  plain  weave  laminate  up  to  200 

hours  using  various  mesh  refinements 


6.5.  Adaptive  time  stepping 

The  adaptive  time  stepping  algorithm  was  validated  by  comparing  results  from  an 
oxidation  analysis  using  adaptive  time  stepping  to  one  that  used  manual  time  stepping. 
Figure  6.19  shows  the  comparison  of  concentration  profiles  at  10,  100,  and  200  hours  for 
a  one-dimensional  oxidation  analysis  of  PMR-15.  The  solid  line  shows  the  result 
obtained  using  adaptive  time  stepping  whereas  the  dashed  line  shows  the  result  from 
manual  time  stepping.  The  adaptive  time  stepping  results  used  a  minimum  time  step 
value  of  10'6  minutes  and  a  adaptive  time  step  tolerance  of  10'3.  The  manual  time 
stepping  results  used  a  uniform  time  step  size  of  0.15  minutes.  While  the  manual  time 
stepping  required  80000  time  steps  to  complete,  the  adaptive  time  stepping  algorithm 
completed  the  analysis  with  only  17,400  time  steps. 
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Figure  6.19  Validation  of  adaptive  time  stepping  algorithm. 


6.6.  Summary 

Various  strategies  were  developed  with  the  goal  of  expediting  oxidation  analysis.  In 
particular,  the  oxidation  analysis  of  textile  composites.  The  complex  architecture  of 
textiles  coupled  with  element  size  constraints  for  oxidation  analysis  can  result  in  highly 
refined  meshes.  Studies  and  validation  were  conducted  to  determine  an  optimal  element 
size  and  time  step  size  for  the  oxidation  of  PMR-15.  Furthermore,  an  adaptive  meshing 
strategy  was  validated  which  allowed  for  an  increase  in  efficiency  by  only  periodically 
solving  the  entire  system  of  equations  in  oxidation  analysis.  A  decoupled  subdomain 
strategy  developed  to  expedite  composite  analysis  was  validated  using  two-dimensional 
and  three-dimensional  models.  Despite  the  significant  simplifying  assumptions,  the 
decoupled  subdomain  strategy  predicts  the  results  of  a  three-dimensional  analysis  very 
well.  The  use  of  exact  and  approximate  boundary  conditions  were  also  examined  for 
three-dimensional  oxidation  analysis  of  a  plain  weave.  It  appears  that  approximate 
boundary  conditions  give  sufficiency  accuracy  and  can  significantly  reduce  the  size  the 
mesh  used  in  oxidation  analysis.  Finally,  an  adaptive  time  stepping  algorithm  was 
implemented  and  validated,  which  showed  very  good  agreement  with  results  obtained  by 
manual  specification  of  an  optimal  time  step. 
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7.  Homogenized  Oxidation  Properties 

When  simulating  oxidation  in  a  configuration  made  of  composites,  it  is  not  practical  to 
discretely  model  all  the  matrix  and  fibers  in  the  composite  because  of  modeling  and 
computational  challenges.  The  same  situation  is  true  in  the  case  of  textile  composites. 
Similar  to  what  is  done  in  order  to  perform  structural  analysis  of  textile  composites, 
homogenized  properties  are  used  to  avoid  modeling  a  microstructural  scale  thereby 
making  the  analysis  tractable.  The  necessity  is  even  more  severe  when  simulating 
oxidation  because  the  finite  element  formulation  requires  very  refined  meshes. 
Therefore,  even  discretely  modeling  the  tow  architecture  scale  in  a  single  unit  cell 
creates  a  very  large  model.  Strategies  for  determining  homogenized  oxidation  properties 
for  unidirectional  laminates  or  tows  are  described  in  Section  3.  In  this  section,  the 
homogenization  strategies  are  validated  using  various  configurations  so  that  they  can  be 
reliably  used  later  to  model  oxidation  in  textile  composites.  The  next  section  describes 
the  material  properties  and  the  configurations  analyzed.  This  is  followed  by  the  results  of 
the  analyses  and  a  discussion  of  the  accuracy  of  the  homogenized  properties. 

7.1.  Material  properties  and  configurations 

Three  configurations  were  analyzed  to  determine  the  accuracy  of  the  homogenized 
oxidation  properties.  The  expressions  described  in  Section  3  were  used  to  determine 
homogenized  oxidation  properties  for  tows  with  a  fiber  fraction  of  28.49%  and  50%. 
Table  5.1  specifies  the  material  properties  for  the  neat  PMR-15  resin  and  Table  7.1 
specifies  those  computed  for  the  homogenized  tows. 

For  a  more  accurate  calculation  of  the  diffusivity  at  Vf=50%,  the  actual  value  of  D  in 
Figure  3.4  obtained  from  micromechanics  (which  is  0.3254)  is  used  rather  than 
calculating  the  value  using  the  formula  for  the  curve  fit  (which  is  0.33).  Mesh  refinement 
and  time  step  sizes  were  determined  such  that  the  analyses  were  computationally 
efficient  while  giving  accurate  results  as  described  in  the  previous  section.  The  three 
configurations  that  were  considered  are  described  next. 
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Table  7.1:  Oxidation  material  properties  for  the  homogenized  tows 


Homogenized  tow 
(Vf=28.49%) 

Homogenized  tow 
(Vf=50%) 

Diffusivity 

Dm 

41.71  xlO'6  mm2/min 
60.87  xlO"6  mm2/min 

34.88  xlO'6  mm2/min 
50.90  xlO'6  mm2/min 

2.50  mol/(m3min) 

1.75  mol/(m3min) 

0OX 

0.4186 

0.5935 

c° 

0.564  mol/m3 

0.395  mol/m3 

GC 

0.01-0.0067(1/40) :  t  <  40 

0.0033  :  t  >  40  (t  in  hours) 

m 

ipc  r  pc  i 

1 +pc[  2(1  +  /?C)J 

p 

0.919 

7.1.1.  Configuration  A  (20  fiber) 

This  configuration  is  a  unidirectional  laminate  idealized  as  having  fibers  arranged  in  a 
“square  array”  with  twenty  fibers  in  the  x  direction  and  infinite  dimensions  in  the  y  and  z 
directions.  This  is  illustrated  in  Figure  7.1  showing  a  single  layer  of  fibers  from  an 
infinite  stack  of  such  layers.  Although  the  sketch  shows  a  finite  z  dimension,  the 
configuration  is  actually  infinite  in  the  z  direction.  The  fibers  are  identical  and  have  a 
diameter  of  10  microns.  The  fiber  volume  fraction  of  the  laminate  is  50%.  The  laminate 
is  exposed  to  air  on  both  the  left  and  right  surfaces.  Therefore,  the  composite  begins 
oxidizing  from  the  outer  surface  with  the  oxygen  making  its  way  into  the  interior  of  the 
laminate.  By  taking  into  consideration  the  symmetries,  this  configuration  can  be 
analyzed  by  modeling  just  ten  fibers  in  a  two-dimensional  model  as  shown  in  Figure  7.1. 
The  analysis  region  is  also  shaded  in  the  sketch  of  the  configuration.  The  left  edge  of  the 
model  is  exposed  to  air  whereas  the  right  edge  is  impermeable.  The  ten  fibers  in  the 
matrix  are  modeled  discretely  and  the  results  from  using  this  model  will  provide  the 
reference  solution.  The  fibers  are  modeled  as  voids  since  the  fibers  are  assumed  to  be 
impermeable.  The  same  configuration  is  analyzed  in  two  other  ways  to  test  the  accuracy 
of  the  effective  properties.  One  is  to  model  the  configuration  completely  using 
homogenized  properties  for  the  microstructure.  Since  this  involves  only  one 
homogenized  material  in  a  simple  one  dimensional  geometry,  this  can  be  analyzed  using 
a  ID  finite  element  model.  The  other  way  is  to  use  a  mixed  model  with  three  unit  cells 
on  the  extremes  modeled  discretely  and  the  four  interior  unit  cells  modeled  using 
homogenized  properties.  Figure  7.2  shows  a  schematic  of  this  model.  This  method  will 
test  the  accuracy  of  the  homogenized  properties  in  models  with  heterogeneous  materials. 
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Figure  7.1:  Schematic  and  analysis  region  for  configuration  A  with  the  numbering 

for  each  unit  cell. 


Homogenized  unit  cells 

Figure  7.2:  Mixed  model  for  configuration  A 


7.1.2.  Configuration  B(  (36  fiber) 

This  configuration  is  slightly  more  complex  than  the  previous  one  in  that  the  laminate  is 
infinite  only  in  the  z  direction.  This  can  be  considered  as  a  square  tow  with  36  fibers 
packed  in  a  square  array  as  illustrated  in  Figure  7.3.  Again,  all  the  fibers  have  a  diameter 
of  10  microns  and  are  packed  with  a  fiber  fraction  of  50%.  The  tow  is  exposed  to  air  on 
all  four  lateral  surfaces  and  starts  oxidizing  as  the  oxygen  makes  diffuses  into  the  tow. 
Utilizing  symmetry  conditions,  only  the  shaded  region  in  the  sketch  needs  to  be 
modeled,  as  shown  in  Figure  7.3.  The  left  and  bottom  edges  in  the  analysis  model  are 
exposed  to  air  while  the  right  and  top  edges  are  specified  to  be  impermeable. 
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Figure  7.3:  Schematic  and  analysis  region  for  configuration  B  with  the 
numbering  for  each  unit  cell. 


7.1.3.  Configuration  C  (irregular  fiber  distribution) 

This  configuration  considers  a  slightly  more  realistic  situation  with  an  irregular 
distribution  of  fibers.  A  computer  generated  random  microstructure  was  used  to 
represent  the  microstructure  in  a  tow  (see  Figure  7.4).  The  two  dimensional 
microstructure  assumes  that  the  fibers  run  exactly  parallel  to  each  other  in  the  z- 
direction.  This,  of  course,  is  not  what  happens  in  a  typical  tow  but  this  configuration 
would  be  a  reasonable  precursor  to  modeling  the  much  more  complex  (if  at  all 
achievable)  realistic  microstructure  of  a  tow.  The  configuration  C  has  160  identical 
fibers  with  10  micron  diameter  like  the  previous  configurations  but  with  an  overall  fiber 
fraction  of  28.49%.  The  analysis  region  is  assumed  to  be  a  square  with  a  side  of  210 
microns.  The  configuration  is  assumed  to  be  exposed  to  air  from  the  bottom  edge  and  all 
the  other  surfaces  are  impermeable.  Three  idealizations  are  used  to  model  the 
configuration.  The  first  one  discretely  models  the  random  microstructure  shown  in 
Figure  7.4.  This  idealization  also  brings  to  light  the  computational  challenges  involved 
in  analyzing  the  oxidization  behavior  of  complex  microstructures. 
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Figure  7.4:  Analysis  regions  for  the  different  configuration  C  idealizations. 

The  second  idealization  uses  a  periodic  microstructure.  It  is  not  possible  to  create  a 
perfect  square  region  using  an  array  of  160  square  unit  cells  because  a/T60  is  not  a 
rational  number.  A  close  approximation  was  chosen  using  a  square  with  a  side 
measuring  12.5  square  unit  cells  (or  207.54  microns).  Although  the  height  of  the 
periodic  model  is  a  half  unit  cell  longer  than  the  discrete  model,  this  difference  does  not 
have  any  effect  on  the  oxidation  growth  behavior  for  the  200  hour  simulations  that  are 
analyzed  in  this  work.  Even  after  200  hours  of  oxidation,  the  oxidation  front  in  a  pure 
resin  advances  less  than  100  microns  and  the  un-oxidized  material  on  the  other  side  of 
the  front  has  insignificant  influence  on  the  oxidation  growth  up  to  that  point.  For  a 
periodic  idealization,  it  is  possible  to  analyze  just  a  fraction  of  the  configuration  by 
taking  advantage  of  symmetry.  To  analyze  the  idealized  periodic  configuration,  a  model 
with  a  width  of  a  half  unit  cell  and  a  height  of  12.5  unit  cells  was  chosen.  The  third 
idealization  uses  a  homogenized  material  to  model  the  configuration.  The  simple 
boundary  conditions  and  the  single  homogenized  material  in  the  idealization  allow  the 
third  configuration  to  be  modeled  using  ID  elements. 


7.2.  Results  and  discussion 

Before  determining  the  accuracy  of  the  effective  oxidation  material  properties,  diffusion 
analyses  were  run  for  both  configurations  A  and  B  to  validate  the  accuracy  of  the 
effective  diffusivities.  Oxidation  analyses  were  conducted  for  all  three  configurations 
(A,  B  and  C).  The  results  from  the  diffusion  modeling  are  discussed  first  followed  by 
that  of  the  oxidation  modeling. 
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7.2.1.  Diffusion  modeling 

The  diffusion  behavior  was  simulated  using  the  un-oxidized  PMR-15  resin  diffusivity  to 
model  the  material  in  the  discrete  models,  which  is  5 3. 6x1  O'6  mm2/min.  For  the  models 
that  used  homogenized  materials,  the  corresponding  effective  diffusivity  of  the  un¬ 
oxidized  resin  was  used,  which  for  the  case  of  a  50%  fiber  fraction  tow  is  34.88  xlO'6 
mm2/min. 

Diffusion  analyses  were  conducted  on  all  three  models  representing  configuration  A: 
discrete  model  which  serves  as  the  reference  solution,  a  fully  homogeneous  model  and  a 
mixed  model  as  shown  in  Figures  7.1  and  7.2.  The  two-dimensional  models  that 
generated  the  results  shown  in  Figures  7.5,  7.6  and  7.7  used  meshes  with  a  maximum 
element  length  of  1.41xl0'3mm  and  a  time  step  size  of  0.15  minutes.  Figure  7.5  shows 
the  concentration  profiles  in  the  discrete  and  mixed  model  at  5  hrs.  It  is  seen  that  the 
concentration  profiles  are  almost  exactly  the  same  in  the  first  three  unit  cells  on  the  left 
which  is  modeled  discretely  in  both  the  discrete  and  mixed  models.  This  shows  that  the 
effective  properties  in  the  homogenized  region  did  not  cause  an  adverse  effect  on  the 
concentration  profile  in  the  discretely  modeled  region.  The  homogenized  material  has  a 
different  saturation  concentration  compared  to  the  neat  PMR-15  material  as  shown  in 
Table  7.1.  The  normalized  concentration  distribution  will  be  continuous  throughout  the 
model  based  on  the  finite  element  formulation  as  explained  in  Section  3.  On  the  other 
hand,  since  the  homogenized  material  has  a  different  saturation  concentration  as 
compared  to  the  neat  resin,  the  concentration  distribution  will  have  a  discontinuity  at  the 
interface  between  the  discrete  and  homogenized  region.  It  is  not  convenient  to  make 
reasonable  comparison  between  the  two  models  by  just  looking  at  the  concentration 
profiles.  When  the  models  compared  have  homogenized  properties,  it  is  perhaps  more 
reasonable  to  compare  volume  averaged  concentrations. 
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Figure  7.5:  Concentration  profiles  in  discrete  and  mixed  models  for  configuration 

A  under  diffusion  at  5  hours. 


Figure  7.6  shows  the  volume  averaged  concentration  for  the  entire  model  as  it  grows 
over  time.  It  shows  that  the  curves  from  the  homogeneous  and  mixed  model  fall  exactly 
on  top  on  the  curve  from  the  discrete  model.  To  take  a  closer  look  at  the  results,  the 
concentrations  were  averaged  over  each  of  the  10  unit  cells  in  the  configuration.  The  unit 
cells  are  numbered  as  shown  in  Figure  7.1.  Figure  7.7  shows  the  average  concentrations 


84 


in  each  unit  cell  at  three  snapshots  in  the  simulation:  15  mins,  150  mins  and  375  mins 
into  the  simulation.  For  each  snapshot,  the  average  concentrations  from  the  three 
different  models  are  shown.  The  first  three  columns  for  each  unit  cell  denote  the  average 
concentrations  for  t=15  mins,  the  next  three  for  t=150  mins  and  the  last  three  columns 
for  t=375  mins.  For  each  set  of  three  columns,  the  first  one  denotes  the  discrete  model, 
the  second  denotes  the  homogeneous  model  and  the  last  one  denotes  the  mixed  model. 
The  results  show  that  both  the  models  that  use  effective  properties  agree  very  well  with 
the  discrete  model. 

Configuration  B  was  analyzed  for  diffusion  using  both  the  discrete  model  and  the  fully 
homogenized  model.  The  two-dimensional  models  that  generated  the  results  shown  in 
Figures  7.8  and  7.9  used  meshes  with  a  maximum  element  length  of  1.41xl0'3mm  and  a 
time  step  size  of  0.15  minutes.  Figure  7.8  shows  the  average  concentration  in  the  entire 
model  as  it  grows  over  time.  The  two  models  agree  closely.  The  homogeneous  model 
under  predicts  the  average  concentration  in  the  beginning  of  the  simulation  and  the 
difference  reduces  as  time  progresses.  It  is  not  surprising  that  the  error  reduces  as  the 
simulation  progresses  because  both  models  are  approaching  the  same  steady  state 
condition.  Similar  to  the  previous  configuration,  the  average  concentration  was 
determined  for  each  of  the  nine  unit  cells  at  two  different  times  through  the  simulation. 
The  unit  cells  are  numbered  as  shown  in  Figure  7.3.  Figure  7.9  shows  the  average 
concentration  from  the  discrete  and  homogeneous  in  each  unit  cell  at  t=1.5  and  15 
minutes.  The  results  in  Figure  7.9  repeat  the  trend  from  Figure  7.8  in  that  the 
homogeneous  model  under  predicts  the  concentrations  and  the  predictions  become  closer 
in  agreement  as  the  simulation  progresses.  In  this  configuration,  the  concentration  profile 
is  more  complicated  than  the  earlier  one  because  the  oxygen  is  flowing  in  from  two 
directions.  This  kind  of  complex  loading  appears  to  have  an  effect  on  how  well  the 
diffusion  is  modeled  by  the  effective  properties.  Although  it  does  not  do  a  perfect  job  in 
simulating  the  oxygen  flow  in  the  beginning,  the  accuracy  increases  very  quickly  as  the 
simulation  progresses. 
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Figure  7.6:  Variation  of  average  concentration  in  configuration  A  with  time 
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Figure  7.7:  Variation  of  average  concentration  in  each  unit  cell  in  configuration  A 

at  different  times  (in  minutes) 
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Figure  7.9:  Variation  of  average  concentration  in  each  unit  cell  in  configuration 

B  at  different  times  (in  minutes) 
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7.2.2.  Oxidation  modeling 

This  section  discusses  the  results  from  the  oxidation  simulation  of  configurations  A,  B 
and  C.  In  the  oxidation  analysis,  there  are  primarily  two  types  of  data  that  are  of  interest 
-  the  concentration  and  the  oxidation  state.  For  configuration  A  and  B,  the  concentration 
distribution  will  be  discussed  first  followed  by  the  oxidation  state.  For  configuration  C, 
only  the  oxidation  state  results  are  presented. 

Figure  7.10  shows  the  concentration  profiles  in  the  discrete  and  mixed  model  at  200  hrs. 
It  is  seen  that  the  concentration  profiles  are  almost  exactly  the  same  in  the  first  three  unit 
cells  on  the  left  which  are  modeled  discretely  in  both  the  discrete  and  mixed  models. 
This  shows  that  the  effective  properties  in  the  homogenized  region  did  not  cause  an 
adverse  effect  on  the  concentration  profile  in  the  discretely  modeled  region.  There  is,  as 
expected,  a  discontinuity  in  the  concentration  at  the  interface  between  the  discrete  and 
homogenized  region  just  as  seen  in  the  diffusion  analysis  of  configuration  A. 
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Figure  7.10:  Concentration  profiles  in  discrete  and  mixed  models  for 
configuration  A  at  200  hours. 


Just  as  the  results  for  the  diffusion  analysis  were  presented,  Figure  7.11  shows  the 
average  concentration  growth  in  the  model  over  time.  The  plot  shows  that  the  discrete 
model  appears  to  have  spurts  of  increase  in  the  average  concentration.  This  can  be 
explained  by  the  fact  that  the  discrete  model  has  fibers  that  are  impermeable  and  do  not 
oxidize.  When  the  oxygen  diffuses  from  the  left  end,  the  cross  sectional  area  of  the 
polymer  material  through  which  it  can  diffuse  varies.  The  area  decreases  to  a  minimum 
where  the  fiber  takes  up  the  most  space  in  the  cross  section  (indicated  by  A  in  Figure 
7.10)  and  increases  to  a  maximum  when  there  is  no  fiber  in  the  cross  section  (indicated 
by  B).  Therefore,  when  the  oxygen  is  diffusing  through  the  constricted  regions,  the 
process  slows  down  and  this  effect  shows  up  in  the  concentration  growth.  When  the 
oxidation  front  (or  the  ‘moving  barrier’  as  described  in  the  previous  section)  passes  the 
constricted  pathways,  the  process  speeds  up  for  a  while  till  the  next  constricted  pathway 
comes  along.  The  homogeneous  model  has  no  such  spurts  in  the  growth  of  the  average 
concentration  because  the  model  assumes  that  it  is  all  one  homogeneous  material 
without  any  impermeable  fibers.  On  the  other  hand,  the  mixed  model  does  show  the 
spurts  in  concentration  growth  in  the  first  part  of  the  curve  because  the  mixed  model 
does  have  three  discrete  fibers  modeled  in  the  left  end  where  the  oxygen  is  entering  the 
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material.  Although  there  are  these  oscillations  when  the  fibers  are  discretely  modeled,  it 
is  clearly  seen  that  the  models  with  the  effective  properties  do  follow  the  same  general 
trend  and  appears  to  follow  the  mean  line  of  the  oscillating  curves. 

Figure  7.12  shows  the  difference  in  the  average  concentration  from  the  discrete  and 
homogeneous  in  each  unit  cell  at  t=150  minutes  and  t=200  hours.  While  the 
homogeneous  models  always  under  predicted  the  average  concentrations  in  the  diffusion 
only  analysis  (see  Figures  7.7  and  7.9),  no  such  correlation  was  seen  in  the  oxidation 
analysis.  This  can  be  attributed  to  the  same  reason  for  seeing  surges  in  the  growth  of  the 
average  concentration.  As  seen  in  Figure  7.11,  depending  on  the  simulation  time,  the 
homogeneous  model  fluctuates  between  under  predicting  and  over  predicting  the 
average  concentration.  This  same  effect  is  what  is  seen  in  Figure  7.12. 

The  Zone  I  and  II  thicknesses  are  measured  for  all  the  three  models  for  configuration  A: 
the  fully  discrete  model,  fully  homogenized  model  and  the  mixed  model.  The  zone 
thicknesses  for  the  discrete  and  mixed  models  are  assumed  to  be  the  thicknesses  along 
the  top  or  bottom  edges  of  the  model.  Note  that  the  model  is  symmetric  about  the 
horizontal  mid-axis  and  therefore  the  oxidation  layer  growth  will  be  symmetric  about 
that  line.  Figure  7.13  shows  growth  of  the  oxidation  layer  (Zone  I  +  II)  for  the  three 
models  as  well  as  the  variation  of  the  active  zone  layer  (Zone  II).  It  can  be  seen  that  the 
effective  properties  do  a  good  job  in  predicting  the  growth  in  both  the  homogenous  and 
mixed  model.  The  Zone  II  thickness  is  also  found  to  be  predicted  fairly  well  considering 
that  the  Zone  II  thickness  according  to  the  discrete  model  appears  to  fluctuate  erratically. 
Figure  7.14  shows  the  evolution  of  the  oxidation  layer  in  the  discrete  and  mixed  models. 
The  three  zones,  Zone  I  (Fully  oxidized),  Zone  II  (Oxidizing)  and  Zone  III  (Un- 
Oxidized)  are  color-coded  by  grey,  red  and  blue  respectively.  Figure  7.14  shows  the 
state  of  oxidation  at  six  different  times  during  the  simulation,  t=  150  mins,  25  hrs,  50 
hrs,  100  hrs,  150  hrs  and  200  hrs.  The  snapshots  in  Figure  7.14  shows  what  has  already 
been  conveyed  by  Figure  7.13  in  that  the  effective  properties  are  able  to  simulate  the 
oxidation  layer  growth  fairly  well  for  configuration  A. 
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Average  Concentration  in  Model  (mol/n) 


Figure  7.1 1:  Variation  of  average  concentration  in  configuration  A  with  time  under 

oxidation 


1  23456789  10 

Unit  Cells 


Figure  7.12:  Variation  of  average  concentration  in  each  unit  cell  in  configuration 
A  at  different  times  under  oxidation. 
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Next,  the  oxidation  analysis  was  performed  for  configuration  B  using  two  models  -  a 
discrete  model  and  a  homogeneous  model.  Figure  7.15  shows  the  concentration 
distribution  in  the  two  models  at  time,  t=150  mins.  The  homogeneous  model  has  the 
fibers  drawn  in  light  gray  in  order  to  expedite  comparisons  with  the  discrete  model. 
While  the  concentrations  contours  do  not  exactly  match,  the  contours  in  between  the 
fibers  do  in  some  sense  resemble  corresponding  contours  in  the  homogeneous  model. 
Figure  7.16  shows  the  growth  of  average  concentration  in  the  two  models  as  simulation 
progresses.  Similar  to  the  concentration  growth  in  discrete  model  for  configuration  A,  it 
can  be  seen  that  the  concentration  growth  for  configuration  B  follows  a  similar 
oscillating  trend.  The  plot  shows  that  the  homogeneous  model  under-predicts  the 
average  concentration  for  most  of  the  simulation.  It  is  understandable  that  the 
homogeneous  model  is  not  able  to  reproduce  the  wavy  nature  of  the  concentration 
growth  in  Figure  7.16  because  that  effect  is  caused  by  the  microstructure.  The 
homogenous  model  on  the  other  hand,  does  not  have  such  a  microstructure. 


Figure  7.13:  Oxidation  layer  growth  (Zone  I+II  and  Zone  I)  for  configuration  A 
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Discrete  Model 


Mixed  Model 


Figure  7.14:  Evolution  of  oxidation  layer  in  discrete  and  mixed  model  for 

configuration  A 
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Figure  7.15:  Concentration  profiles  in  discrete  and  mixed  model  at  t=150  mins 

for  configuration  A 
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Figure  7.16:  Variation  of  average  concentration  in  configuration  B  with  time  under 

oxidation 


The  oxidation  state  from  the  two  models  is  compared  next.  Figure  7.17  shows  the 
oxidation  state  at  different  times  in  the  simulation.  As  shown  in  Figure  7.15,  the 
homogeneous  model  has  the  fibers  drawn  in  light  gray  in  order  to  expedite  comparisons 
with  the  discrete  model.  The  oxidation  states  for  t=15  mins,  2.5  hrs,  5  hrs  and  10  hrs  are 
shown.  The  figures  show  that  the  homogeneous  model  is  able  to  predict  the  oxidation 
layer  growth  fairly  accurately. 

The  discrete  models  that  have  been  considered  so  far  model  only  9  or  10  fibers  but  the 
discrete  model  for  configuration  C  models  160  fibers.  This  makes  it  a  considerably 
larger  model  and  more  time-consuming  compared  to  the  previous  models.  Figure  7.18 
shows  the  contour  plots  of  the  oxidation  state  in  the  discrete  and  periodic  models  after 
undergoing  75  hours  of  oxidation.  As  expected  the  oxidation  profile  is  irregular  for  the 
model  with  the  discrete  microstructure.  Nevertheless,  the  variation  in  thickness  and 
location  of  the  oxidation  layer  from  the  exposed  edge  varies  no  more  than  5%  across  the 
width.  It  is  interesting  to  note  that  the  periodic  model  predicts  relatively  the  same 
amount  of  oxidation  growth  as  the  random  model  which  indicates  that  for  this  fiber 
volume  fraction  and  distribution  of  fibers,  the  oxidation  growth  can  be  idealized  by 
using  a  periodic  array.  While  Figure  7.18  shows  the  oxidation  state  distribution  at  a 
single  snapshot  from  the  entire  simulation,  Figure  7.19  gives  a  sense  of  how  the 
oxidation  state  evolved  during  the  entire  simulation.  Figure  7.19  gives  the  oxidation 
layer  growth  over  time  for  the  random,  periodic  and  homogenized  idealizations.  It  shows 
the  oxidation  layer  growth  along  the  two  edges  (right  and  left)  of  the  discrete  model. 
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While  the  two  curves  do  not  fall  right  on  top  of  each  other,  they  are  very  close.  The 
curve  from  the  periodic  model  is  very  close  to  the  curves  from  the  discrete  model  and 
follows  the  same  trend  but  slightly  under  predicts  the  oxidation  growth.  The 
homogeneous  model  also  follows  the  same  trend  but  under  predicts  the  growth  even 
further.  To  make  an  easy  comparison  with  the  behavior  if  there  were  no  fibers  at  all  (i.e. 
pure  resin),  the  oxidation  growth  curve  from  pure  resin  oxidation  analysis  (using  a  ID 
model)  is  also  included.  This  shows  that  the  pure  resin  oxidizes  slightly  faster  than  when 
there  are  fibers  in  the  resin  which  is  expected  since  the  fibers  are  assumed  to  be 
impermeable  and  do  not  oxidize. 


t=  l  5  mins 

Discrete  Homogeneous 


Discrete  Homogeneous 

t=5  hrs 


t=2.5  hrs 

Discrete  Homogeneous 


Discrete  Homogeneous 

t=10  hrs 


Figure  7.17:  Evolution  of  oxidation  layer  in  discrete  and  homogeneous  models  for 

configuration  B 
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Oxidation  Layer  Thickness  (mm) 


Random  Periodic 

microstructure  microstructure 


Figure  7.18:  Oxidation  state  profiles  in  discrete  and  periodic  model  at  t=75 

hours  for  configuration  C 
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A  typical  homogenization  process  in  structural  mechanics  results  in  being  able  to  use  a 
less  refined  homogenized  model  to  replace  the  actual  microstructure.  This  was  generally 
found  to  be  the  case  for  the  oxidation  analysis  as  well.  As  discussed  in  the  previous 
section,  it  is  possible  optimize  the  mesh  parameters  and  use  larger  elements  and  time 
steps.  In  some  configurations,  including  the  ones  described  in  this  section,  the  element 
size  is  restricted  by  the  need  model  the  discrete  geometry  accurately.  This  restriction  is 
greatly  reduced  when  modeling  homogenized  regions.  It  is  also  easier  to  generate  the 
models  when  a  complex  microstructure  can  be  replaced  by  a  simpler  homogenized 
geometry.  Another  advantage  is  that  sometimes  a  two  dimensional  model  can  be 
replaced  with  a  one  dimensional  model  that  is  much  less  computationally  intensive.  This 
was  made  use  of  when  analyzing  configuration  A  with  a  fully  homogeneous  model.  The 
goal  of  this  work  is  to  eventually  be  able  to  correlate  the  effect  of  the  oxidation  on  the 
mechanical  response  and  damage  progression  of  the  composite.  Keeping  this  in  mind, 
tracking  the  regions  of  oxidized  or  oxidizing  material  in  the  composite  is  what  would  be 
considered  to  impact  the  mechanical  response.  For  the  material  system  being  considered 
in  this  work,  the  thickness  of  the  active  zone  is  found  to  be  practically  constant  and  the 
variation  of  oxidation  state  within  this  zone  can  be  inconsequential  in  this  regard,  but 
this  need  not  be  the  case  for  other  composite  systems.  Further  work  needs  to  be 
performed  in  order  to  determine  if  some  accuracy  of  the  oxidation  state  variation  in  the 
active  zone  can  be  given  up  in  exchange  for  better  computational  efficiency  as  long  as 
the  oxidation  layer  thicknesses  are  predicted  with  reasonable  accuracy. 

7.3.  Summary 

To  simulate  the  oxidation  of  the  textile  composite,  it  is  important  to  be  able  to  use 
homogenized  oxidation  properties  for  the  tow  because  it  is  practically  impossible  to 
discretely  model  all  the  fibers  in  a  composite.  Effective  oxidation  material  properties 
were  calculated  for  a  unidirectional  laminate/tow  using  the  expressions  described  in 
Section  3.  Three  configurations  were  analyzed  to  test  the  accuracy  of  the  effective 
oxidation  properties.  The  fibers  were  assumed  to  be  impermeable  and  do  not  oxidize.  All 
the  configurations  had  10  micron  diameter  circular  fibers.  Two  of  the  configurations  had 
the  fibers  in  a  square  array  packing  with  50%  fiber  fraction  whereas  the  third 
configuration  had  random  microstructure  with  an  overall  fiber  fraction  of  28.5%.  The 
configurations  were  discretely  modeled  to  provide  a  reference  solution.  Idealizations 
with  fully  homogenized  materials  as  well  as  mixed  idealizations  (both  discrete  and 
homogenized  regions)  were  used  to  determine  the  accuracy  of  the  effective  properties. 
The  concentration  of  oxygen  in  the  model  as  well  as  the  oxidation  state  of  the  materials 
in  the  composite  was  compared  to  the  reference  model.  It  was  seen  that  the  effective 
oxidation  properties  performed  reasonably  well  for  both  configurations  and  were  able  to 
simulate  the  oxidation  layer  growth. 
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8.  Oxidation  Analysis  of  Textile  Composites 

One  of  the  primary  goals  of  this  work  is  to  study  the  effect  of  oxidation  on  damage 
progression  in  textile  composites.  A  precursor  to  the  complete  damage  progression  is  the 
oxidation  analysis  of  the  textile  composite.  The  information  from  the  oxidation  analysis 
will  be  used  to  degrade  the  mechanical  properties  of  the  textile  composite  in  the  damage 
progression  model.  In  this  work,  the  mechanical  damage  is  assumed  to  not  have  an  effect 
on  the  oxidation  behavior.  Therefore,  the  oxidation  model  and  the  damage  progression 
model  are  only  coupled  in  one  direction,  where  the  oxidation  behavior  has  an  effect  on 
the  damage  model  and  not  the  other  way  round.  As  discussed  in  the  previous  sections, 
simulating  the  oxidation  behavior  is  a  computationally  intensive  task.  At  the  time  the 
results  in  this  section  were  obtained,  the  efficiency  of  the  software  did  not  allow  for  3D 
results  to  be  obtained.  These  results  were  obtained  using  the  decoupled  subdomain 
strategy  (DSS)  described  in  Section  5  and  validated  in  Section  6.  Validation  showed  that 
DSS  agrees  quite  well  with  the  3D  results.  Thus,  the  results  presented  in  this  section 
should  be  viewed  as  having  adequate  accuracy  despite  the  simplifications  imposed  by 
DSS.  This  section  also  describes  the  information  that  is  generated  from  DSS,  which  can 
then  be  used  in  the  coupled  damage  progression  model  to  predict  the  mechanical 
behavior  of  the  composite  under  oxidation. 


8.1.  Plain  weave  analyzed  using  decoupled  subdomain  strategy 

Now  that  DSS  has  been  validated  in  Section  6,  it  can  be  used  to  simulate  the  oxidation 
behavior  in  textile  composites  and  eventually  coupled  with  the  damage  progression 
analysis  to  predict  the  mechanical  behavior  under  oxidation.  The  configurations  that  are 
examined  in  this  work  are  plain  weave  laminates  exposed  to  oxygen  at  the  top  and 
bottom  surfaces  and  the  oxidation  behavior  is  simulated  only  for  200  hours.  As 
discussed  later  in  this  section,  after  200  hours  of  oxidation  of  the  laminates  with  the 
material  system  that  is  considered  in  this  work,  the  oxidation  layer  thickness  does  not 
exceed  more  than  the  thickness  of  a  single  ply.  Therefore,  based  on  the  oxidation  model 
implemented  in  this  work,  the  growth  of  the  oxidation  layer  would  be  the  same 
regardless  of  whether  it  is  a  2-ply  laminate  or  if  it  has  more  than  2  plies.  On  the  other 
hand,  although  the  oxidation  layer  growth  is  the  same,  the  number  of  plies  could  have  an 
impact  on  the  mechanical  behavior  and  this  is  discussed  in  the  next  section. 

DSS  was  used  to  simulate  the  oxidation  behavior  in  a  symmetric  two-ply  graphite/PMR- 
15  plain  weave  laminate.  Both  the  top  and  bottom  surfaces  are  exposed  to  oxygen.  The 
composite  is  chosen  to  have  a  waviness  ratio  of  1/3.  A  full  unit  cell  of  the  configuration 
is  shown  in  Figure  8.1(a).  By  exploiting  symmetry,  it  is  possible  to  analyze  the 
configuration  using  only  l/8th  of  the  full  unit  cell  as  shown  in  Figure  8.1(b)  with  a 
transparent  matrix.  The  decoupled  subdomain  strategy  is  used  on  the  reduced  domain, 
which  is  automatically  subdivided  into  an  array  of  64  ID  models  as  described  in  the 
previous  section.  Since  the  effects  of  the  undulation  are  assumed  to  be  insignificant, 
both  the  warp  and  fill  tows  have  the  same  oxidation  material  properties  in  the  through- 
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thickness  direction.  This  means  that  the  four  quadrants  in  Figure  8.1(b)  are  essentially 
identical  after  the  appropriate  rotations  about  the  z-axis  ,  and  therefore  the  results  from 
the  corresponding  ID  models  in  the  different  quadrants  will  be  the  same.  Figure  8.1(c) 
shows  one  of  the  quadrants  from  Figure  8.1(b)  (labeled  1)  which  is  l/32nd  of  the  full  unit 
cell.  The  labeling  in  Figure  8.1(c)  shows  the  unique  domains  once  the  decoupled 
subdomain  strategy  is  applied.  Since  through  thickness  properties  of  the  tows  are  the 
same  when  the  undulation  effect  on  material  orientation  is  disregarded,  a  plane  of 
symmetry  exists  as  shown  in  Figures  8.1(c)  and  8.1(d).  Figure  8.1(d)  provides  a  view 
rotated  by  180°  about  the  z-axis  to  verify  the  available  symmetry  within  the  l/32nd  unit 
cell  when  employing  the  decoupled  subdomain  strategy. 


(c)l/32nd  unit  cell  with  unique  domains  labeled  (d)  Alternate  view  of  l/32nd  unit  cell  to  display 

symmetries 

Figure  8.1:  Configuration  and  analysis  domains  for  simulating  oxidation  in  plain 

weave  composite 


Therefore,  the  only  unique  results  from  the  analysis  are  those  from  the  10  domains 
numbered  in  Figure  8.1(d).  Figure  8.2  gives  the  predicted  oxidation  layer  growth  for  the 
10  domains.  It  shows  that  there  is  considerable  variation  in  the  oxidation  layer  growth 
behavior  of  the  10  domains.  At  the  end  of  200  hours  of  oxidation,  the  thickest  layer  is 
0.11  mm  (in  Domain  9)  which  is  only  slightly  larger  than  half  the  thickness  of  a  single 
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ply.  Figure  8.1(c)  shows  that  domain  10  has  the  largest  amount  of  matrix  with  a  very 
small  region  of  tow  in  the  middle  whereas  domain  1  has  the  largest  amount  of  tow  with  a 
very  small  region  of  matrix  at  the  two  ends.  Although  domain  10  has  the  largest  amount 
of  matrix,  it  is  not  the  domain  that  has  the  thickest  oxidation  layer.  This  is  because,  as 
discussed  earlier,  in  heterogeneous  models  the  oxidation  behavior  is  not  very 
straightforward  and  depends  on  a  number  of  factors  such  as  the  location  of  the  material 
boundaries  and  the  oxidation  properties  of  each  of  the  constituent  materials.  In  each  of 
the  ten  unique  1-D  domains  representing  the  weave  microstructure,  the  material 
boundaries  are  at  a  different  distance  away  from  the  exposed  surface.  This  results  in  a 
varied  oxidation  behavior  from  the  1-D  models.  Since  domain  10  is  almost  all  resin  with 
a  small  region  of  tow  in  the  middle,  its  oxidation  behavior  would  be  expected  to  be  close 
to  that  of  a  neat  resin.  Similarly,  since  domain  1  is  almost  all  tow  with  small  regions  of 
matrix  at  the  two  ends,  its  oxidation  behavior  would  be  expected  to  be  close  to  that  of  a 
homogenized  tow  model. 

However,  as  explained  earlier  with  the  heterogeneous  configuration,  the  behavior  is  not 
always  close  to  that  of  the  corresponding  homogeneous  model.  Figure  8.3  shows  the 
layer  growth  for  domains  1  and  10  as  well  as  for  a  neat  resin  model  and  a  homogenized 
tow  model.  It  shows  that  domain  10  follows  the  same  behavior  as  a  neat  resin  model  but 
once  the  oxidation  front  reaches  the  tow  material,  domain  10  has  a  slightly  faster 
oxidation  layer  growth.  For  domain  1,  the  oxidation  layer  is  only  slightly  thicker  than 
that  in  an  all  tow  model.  Overall,  the  analysis  shows  that  the  oxidation  front  does  not 
advance  uniformly  throughout  the  composite.  At  the  end  of  200  hours  of  oxidation, 
domain  1  has  the  smallest  oxidation  layer  with  a  thickness  of  84.5  microns  and  domain  9 
has  the  largest  oxidation  layer  with  a  thickness  of  110  microns.  That  is  a  range  of  over 
25  microns,  which  is  over  30%  of  the  domain  1  layer  thickness.  Therefore,  it  can  be  said 
that  the  tow  architecture  plays  a  significant  role  in  the  variation. 
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oxidation  layer  thickness  (mm)  oxidation  layer  thickness 


Figure  8.1:  Oxidation  layer  growth  in  the  10  unique  domains 


Figure  8.2:  Comparison  of  oxidation  layer  growth  in  domains  1  and  10  with  that 
of  a  neat  resin  model  and  homogenized  tow  model 
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8.2.  Storage  of  oxidation  behavior  data  from  decoupled  subdomain 
strategy 

The  oxidation  behavior  of  the  laminate  is  eventually  used  in  the  coupled  damage 
progression  model  in  order  to  predict  the  mechanical  behavior  under  oxidation.  In  order 
to  do  this,  the  results  from  the  oxidation  analysis  need  to  be  passed  on  to  the  damage 
progression  model.  The  results  consist  of  the  distribution  of  the  oxidation  level  property 
in  the  laminate  at  different  time  steps  in  the  simulation.  The  oxidation  level  data  at  the 
different  time  steps  is  needed  by  the  damage  progression  model  in  order  to  degrade  the 
mechanical  properties  of  the  composites  based  on  how  much  of  the  material  has 
oxidized.  The  value  of  the  oxidation  level  at  each  integration  point  in  the  all  the  elements 
of  the  finite  element  model  is  kept  track  of  in  the  memory  and  can  be  written  to  a  file, 
similar  to  how  the  stress  distribution  in  a  model  can  be  written  to  a  file.  If  the  oxidation 
model  and  the  damage  progression  model  used  the  same  discretization  for  the  analysis 
domain,  i.e.  the  same  finite  element  mesh,  then  the  information  transfer  is 
straightforward.  The  oxidation  level  distribution  file  can  be  read  in  by  the  damage 
progression  model  and  all  the  oxidation  level  information  would  be  available  for 
performing  the  mechanical  property  degradation. 

However,  when  a  hybrid  model  is  used  for  the  oxidation  analysis,  the  information 
transfer  to  the  damage  progression  model  is  not  so  straightforward.  In  the  hybrid  model, 
each  individual  1-D  domain  is  an  approximation  of  the  actual  3-D  region  that  it 
represents  in  the  laminate.  Due  to  this  reason,  the  oxidation  level  value  distribution  in 
the  ID  model  is  not  an  exact  representation  of  what  the  distribution  would  be  if  the 
actual  3-D  domain  was  analyzed.  For  example,  Figure  8.4  shows  a  3D  domain  and  its 
equivalent  ID  domain.  Point  A  in  the  3D  domain  would  be  the  geometrically  equivalent 
point  to  Point  B  in  the  equivalent  ID  domain,  but  note  that  the  two  points  are  located  in 
different  material  regions  of  the  models.  Point  A  is  located  in  the  matrix  region  where  as 
Point  B  is  located  in  the  tow  region.  However,  because  of  the  characteristic  oxidation 
behavior,  the  mismatch  in  the  geometry  is  only  an  issue  when  the  oxidation  front  is  in 
the  vicinity  of  the  material  boundaries.  Even  when  the  oxidation  front  is  near  the 
material  boundary,  it  is  seen  that  errors  due  to  this  mismatch  are  not  significant  because 
the  rotation  angles  of  the  tow  in  the  laminates  are  not  large  enough. 

Another  issue  has  to  do  with  the  amount  of  information  that  has  to  be  transferred  from 
the  oxidation  model  to  the  damage  progression  model.  For  example,  using  DSS  on  the 
3D  domain  shown  in  Figure  8.1(b)  would  result  in  64  ID  domains.  Assuming  that  the 
oxidation  level  distribution  is  post-processed  and  outputted  by  each  ID  model  at  88 
different  time  steps  throughout  the  oxidation  simulation,  this  would  result  in  the  creation 
of  88  x  64  =  5632  data  sets.  During  the  coupled  simulation,  at  each  of  the  88  time  steps, 
64  different  files  need  to  be  opened  and  the  oxidation  level  information  of  each  element 
in  all  the  ID  models  need  to  be  read  in.  A  strategy  was  sought  that  could  reduce  the 
number  of  file  I/O  operations  as  well  as  the  amount  of  data  that  had  transferred  during 
the  coupled  analysis  while  maintaining  reasonable  accuracy.  The  approximations  that 
were  made  in  the  developed  strategy  and  a  description  of  the  data  that  is  transferred  from 
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the  hybrid  oxidation  model  to  the  damage  model  are  described  in  the  remainder  of  this 
section. 


Figure  8.3:  3D  domain  and  equivalent  ID  domain  in  hybrid  modeling  strategy 


Note  that  the  value  of  the  oxidation  level  at  a  material  point  can  vary  from  1  to  0  as 
described  in  Section  3.  A  value  of  1  denotes  that  the  material  is  un-oxidized  and  a  value 
of  0  denotes  that  the  material  is  fully  oxidized.  Typically,  a  significant  majority  of  the 
model  is  made  up  of  either  fully  oxidized  or  un-oxidized  material.  A  small  fraction  of 
the  model  has  oxidation  levels  in  the  range  between  1  and  0,  which  ideally  denotes  the 
active  zone,  or  that  the  material  has  started  oxidizing  but  it  is  not  fully  oxidized  yet  as 
shown  in  Figure  3.5.  Therefore,  instead  of  storing  the  oxidation  level  information  for 
each  element  in  the  1-D  model,  just  the  dimensions  of  the  active  zone  is  stored  to 
represent  the  oxidation  level  profile  for  a  particular  time  step.  In  this  manner,  the 
oxidation  level  profile  in  a  ID  domain  for  all  of  the  88  time  steps  can  be  effectively 
compressed  into  a  single  file  using  only  a  fraction  of  the  information.  When  the  data  is 
read  in  during  the  coupled  simulation,  the  oxidation  level  profile  is  approximated  using  a 
linear  variation  of  the  oxidation  level  within  the  active  zone.  These  approximations  are 
made  based  on  a  few  assumptions.  The  active  zone  is  assumed  to  be  very  small 
compared  to  the  fully  oxidized  and  un-oxidized  region.  Although  the  actual  variation  of 
the  oxidation  level  in  the  active  zone  is  not  linear,  the  linear  variation  assumed  in  this 
model  is  assumed  to  be  reliable  for  the  material  systems  considered  in  this  work.  The 
simple  linear  approximations  employed  here  are  assumed  to  be  reliable  for  the  purposes 
of  predicting  mechanical  behavior  in  the  composites.  Figure  8.5  shows  the  predicted 
oxidation  level  profile  in  a  block  of  neat  resin  at  100  hours.  The  dotted  line  shows  the 
approximated  oxidation  level  profile.  The  location  where  the  approximated  oxidation 
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Position  (mm) 


Figure  8.4:  Oxidation  Level  profile  in  neat  resin  1-D  model  at  100  hours 


level  starts  lowering  from  1.0  is  determined  by  the  thickness  of  the  oxidation  layer  or  in 
other  words  the  dimensions  of  the  oxidation  zones.  The  instructions  to  determine  the 
dimensions  of  the  different  zones  and  the  oxidation  layer  thickness  are  described  in 
Section  3.  The  location  of  the  point  where  the  approximated  oxidation  level  reaches  0  is 
also  similarly  determined  by  the  dimensions  of  the  active  zone  (typically  it  is  the 
location  where  the  predicted  oxidation  level  reaches  0.01).  During  the  initial  stages  of 
oxidation  when  there  is  no  fully  oxidized  material,  the  predicted  oxidation  level  does  not 
drop  all  the  way  to  0.  In  this  case,  the  linear  approximation  is  based  on  the  value  of  the 
predicted  oxidation  level  at  the  location  in  the  model  that  is  exposed  to  oxygen.  This  is 
illustrated  in  Figure  8.6,  which  shows  the  oxidation  level  profile  in  neat  resin  at  1  hour. 
The  predicted  oxidation  level  at  the  exposed  surface  after  1  hour  of  oxidation  is  0.1812 
and  as  shown  in  Figure  8.5,  both  the  predicted  profile  and  the  approximate  profile  have 
the  same  oxidation  level  value  at  the  exposed  end.  When  analyzing  heterogeneous 
models,  the  oxidation  level  profile  is  more  complicated  in  that  the  profile  is  piece-wise 
continuous  with  the  predicted  oxidation  level  continuous  within  a  single  material.  For 
example,  in  the  heterogeneous  configuration  shown  in  Figure  8.4,  the  material  boundary 
is  at  0.06  mm.  Figure  8.7  shows  the  predicted  oxidation  level  profile  for  that 
configuration  at  70  hours.  The  approximated  oxidation 
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Figure  8.5:  Oxidation  Level  profile  in  neat  resin  1-D  model  at  1  hour 


level  profile  is  also  maintained  as  a  piece-wise  oxidation  level  profile  for  each  material 
region.  The  approximated  oxidation  level  value  in  either  material  region  at  the  material 
boundary  is  same  as  the  corresponding  predicted  oxidation  level  value  for  that  location. 

In  order  to  save  the  approximated  oxidation  level  profiles  for  the  required  time  steps  so 
that  it  can  be  used  by  the  damage  progression  model  in  the  coupled  simulation,  the 
oxidation  level  information  for  each  material  region  in  the  ID  domain  is  stored  using 
just  four  values  -  the  beginning  and  end  locations,  and  the  beginning  and  end  oxidation 
level  values.  This  information  is  then  used  in  the  coupled  damage  progression  model  to 
determine  the  oxidation  level  at  each  integration  point  and  degrade  the  mechanical 
properties  based  on  the  constitutive  relations. 
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Figure  8.6:  Oxidation  Level  profile  in  heterogeneous  1-D  model  (see  Figure  8.4) 
at  70  hours 


8.3.  Summary 

Finite  element  analysis  of  the  oxidation  of  a  plain  weave  graphite/PMR-15  textile  composite  was 
performed  using  the  decoupled  subdomain  strategy  (DSS).  This  strategy  was  presented  and 
validated  in  previous  sections.  This  section  presented  the  oxidation  layer  growth  through  the 
thickness  of  a  woven  mat  by  inspecting  the  results  of  individual  ID  domains.  The  section  also 
describes  how  the  oxidation  information  from  the  DSS  model  is  maintained  so  that  it  can 
be  used  by  the  damage  progression  model  for  prediction  of  mechanical  behavior  of  a 
composite  subjected  to  a  high  temperature  oxidizing  environment. 
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9.  Oxidation-Damage  Analysis  of  Textile  Composites 

The  previous  sections  have  laid  the  groundwork  in  order  to  be  able  to  build  the 
framework  so  that  the  effect  of  oxidation  on  the  mechanical  behavior  of  textile 
composites  can  be  predicted.  The  last  three  sections  describe  the  challenges  and 
appropriate  strategies  for  simulating  the  oxidation  behavior  in  textile  composites. 
Section  3  described  the  governing  equations  and  finite  element  formulations  required  for 
the  damage  analysis,  oxidation  analysis  and  the  coupled  oxidation-damage  progression 
model.  This  section  begins  with  a  brief  overview  of  the  different  damage  mechanisms  in 
textile  composites.  This  is  followed  by  a  description  of  how  the  coupled  analysis  model 
was  used  to  predict  the  damage  initiation  and  progression  in  the  textile  composites  in 
oxidizing  environments.  The  configurations  that  will  be  analyzed  will  be  described 
including  the  material  properties  and  the  constitutive  model  that  was  used  to  implement 
the  coupled  analysis  model.  This  will  be  followed  by  the  results  and  discussion  of  the 
analysis  and  the  parametric  studies. 


9.1.  Damage  mechanisms 

Textile  composites  fail  under  different  types  of  loadings  exhibiting  different  types  of 
damage  mechanisms  [48].  One  common  damage  mechanism  is  transverse  cracking  in 
the  matrix  and  tows.  Other  damage  mechanisms  seen  in  the  tows  are  inter-  and  intra-tow 
delamination,  fiber  buckling  and  fiber  breakage  etc.  Resin  pockets  in  the  composite  can 
develop  transverse  matrix  cracks  (transverse  to  the  loading  direction)  independent  of  the 
matrix  cracks  in  the  tows.  Quaresimin  et  al.  [48]  observed  three  main  damage 
mechanisms  in  twill  weave  composites  under  fatigue  loading.  They  are  layer 
delaminations,  transverse  matrix  cracking  and  fiber  failure.  Figure  9.1  shows  the 
micrographs  illustrating  these  damage  mechanisms.  Quaresimin  et  al.  [48]  analyzed  a 
number  of  laminates  with  different  stacking  sequences  and  saw  that  all  three  damage 
mechanisms  were  present  under  different  types  of  fatigue  loadings,  but  the  sequence  of 
appearance  was  different.  It  was  also  seen  that  only  one  predominant  mechanism 
generally  dictated  the  laminate  behavior. 

Figure  9.2  shows  a  schematic  of  the  different  damage  modes  in  a  tow.  The  mode  under 
which  damage  occurs  in  the  material  depends  on  which  material  allowable  is  exceeded. 
The  failure  criteria  that  are  used  in  this  work  are  discussed  in  the  next  section.  The 
damage  modes  in  the  tows  can  be  classified  into  four  types  as  shown  in  Figure  9.2.  The 
“1”  direction  denotes  the  fiber  direction  whereas  the  “2”  and  “3”  directions  are  in-plane 
and  out  of  plane  transverse  directions  respectively.  The  coordinate  axes  defined  by  the 
“1”,  “2”  and  “3”  direction  are  the  principal  coordinate  axes  of  the  tow,  which  is  assumed 
to  be  transversely  orthotropic.  The  finite  element  model  of  the  composite  accounts  for 
the  undulation  of  the  tows  and  therefore  the  rotation  angles  for  the  material  vary 
depending  on  the  location  of  the  material  point  in  the  tow.  As  illustrated  in  Figure  9.2, 
fiber  breakage  occurs  under  failure  mode  1 1  and  this  damage  mode  is  generally  caused 
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by  excessive  crn  stress  in  the  tow.  This  failure  mode  is  what  generally  causes  the 

ultimate  failure  of  the  composite.  Transverse  matrix  cracking  is  generally  one  of  early 
damage  mechanisms  seen  in  the  tows.  This  type  of  damage  mode  is  caused  by  excessive 
<7 21  or  cr12  stress  components  and  classified  as  failure  modes  22  and  12  respectively. 

Failure  mode  33  and  13  can  be  caused  by  either  <r33  or  a13  stress  components  and  can 

result  in  intra-  or  inter-laminar  delaminations.  Figure  9.2  also  shows  the  damage  mode 
23  which  is  caused  by  <r23  stress. 


(a)  Delamination  (b)  Transverse  matrix  crack 


(c)  Fiber  Failure 

Figure  9.1:  Damage  Mechanisms  in  woven  composites  [48] 
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Figure  9.2:  Schematic  of  different  damage  modes  in  the  tow  of  textile 
composites  [49] 


9.2.  Failure  criteria  for  tows  and  matrix 

As  mentioned  in  the  previous  section,  the  condition  for  damage  to  occur  and  more 
specifically,  which  type  of  damage  mechanism  is  in  action,  is  determined  based  on  what 
failure  mode  has  been  triggered.  In  order  to  determine  if  a  failure  mode  has  been 
triggered,  a  suitable  failure  criterion  is  required.  This  section  defines  the  failure  criteria 
employed  in  all  the  models  used  in  this  work. 

Depending  upon  the  property  degradation  scheme  used,  a  material  point  in  the  matrix 
will  be  assumed  to  be  isotropic  or  anisotropic  after  the  damage  has  occurred.  Since  the 
matrix  is  initially  isotropic,  the  global  coordinate  system  and  the  material  coordinate 
system  are  the  same.  On  the  other  hand,  the  principal  coordinate  system  is  not 
necessarily  the  same  as  the  global  coordinate  system.  Moreover,  the  property 
degradation  scheme  used  in  this  work  assumes  that  the  material  will  become  anisotropic 
after  mechanical  damage.  Therefore,  it  was  assumed  that  there  is  no  significant  effect  of 
choosing  the  maximum  stress  criterion  over  the  principal  stress  criterion.  In  this 
particular  work,  the  stress  in  the  global  coordinate  system  was  used  in  the  maximum 
stress  criterion  to  determine  failure  modes  in  the  matrix.  However,  future  enhancements 
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to  the  model  should  provide  the  option  of  choosing  the  maximum  principal  stress 
criterion  if  the  material  is  not  damaged. 

In  the  case  of  the  tow  material,  the  maximum  stress  criterion  for  anisotropic  materials 
was  used,  which  says  that  the  failure  occurs  when  any  of  the  stress  components  in  the 
material  coordinates  system  exceeds  its  corresponding  strength.  The  tows  can  fail  under 
one  or  more  damage  modes  such  as  fiber  breaking  and  transverse  cracking.  The  modes 
strongly  affect  the  mechanical  behavior  of  the  structure.  In  this  work,  the  tows  are 
assumed  to  be  transversely  isotropic  before  any  damage  occurs.  However,  the  tow  in 
general  would  no  longer  be  transversely  isotropic  after  it  has  failed  and  its  mechanical 
properties  have  been  degraded.  But  the  stress  in  the  local  coordinate  system  of  the  tow  is 
continued  to  be  used  for  the  maximum  stress  failure  criterion.  If  cr,  are  the  stress 

components  in  the  material  coordinates  system  of  the  tow  and  Stj  are  the  corresponding 
strength  values,  then  the  failure  criteria  used  in  this  work  can  be  summarized  as  below: 
For  isotropic  matrix: 

Max  stress  criterion 


For  transversely  isotropic  tow: 
Max  stress  criterion 


>1  Material  point  has  failed  in  mode  ij 


<x. 


''/  <  l  Material  point  has  not  failed 


(6.1) 


The  strength  properties  of  the  materials  analyzed  in  this  work  is  listed  in  Section  9.4 
which  defines  the  all  material  properties  used  in  this  work. 


9.3.  Configuration 

The  coupled  analysis  framework  was  used  to  investigate  the  mechanical  behavior  of  a 
plain  weave  Graphite/PMR-15  composite  in  an  oxidizing  environment.  The  waviness 
ratio  of  the  composite  is  assumed  to  be  1/3  and  the  fiber  volume  fraction  in  the  tow  is 
assumed  to  be  55.6%.  The  reason  for  choosing  this  fiber  volume  fraction  was  because  of 
the  availability  of  experimental  engineering  properties  for  this  particular  material  system 
in  the  literature.  The  volume  fraction  of  the  tows  in  the  composite  is  assumed  to  be 
63.6%  and  therefore  the  overall  fiber  fraction  in  the  composite  is  35.36%.  The  laminate 
consists  of  two  plies  and  is  assumed  to  be  symmetrically  stacked  and  infinite  in  the  in- 
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plane  directions.  Figure  9.3  shows  the  full  unit  cell  of  the  configuration.  The  laminate  is 
assumed  to  be  loaded  under  uniaxial  stress  conditions  in  the  x-direction  at  a  temperature 
of  288C.  The  material  properties  used  to  model  the  configuration  will  be  assumed  to  be 
those  at  288  C.  However,  in  this  current  implementation  of  the  coupled  analysis  model, 
the  effects  of  thermal  expansion  and  the  shrinkage  of  resin  under  oxidation  will  be 
ignored.  The  configuration  will  be  loaded  to  a  pre-determined  strain  level  and 
maintained  at  that  level.  The  configuration  is  then  exposed  to  oxygen  from  the  top  and 
bottom  surfaces  while  the  lateral  surfaces  are  assumed  to  be  impermeable.  The  laminate 
will  be  exposed  to  the  oxygen  for  200  hours  at  28  8C.  The  damage  in  the  laminate 
throughout  this  simulation  will  be  tracked  and  the  mechanical  behavior  will  be  recorded. 
The  number  of  plies  in  the  laminate  was  changed  in  a  parametric  study  to  determine  its 
effect  on  the  mechanical  behavior.  Any  changes  from  this  basic  configuration  will  be 
described  as  required  when  discussing  the  results  of  the  parametric  study. 


Figure  9.3:  Two-ply  plain  weave  composite  configuration 


Exploiting  symmetry  conditions  in  this  configuration  allows  reducing  the  analysis 
domain  from  a  full  unit  cell  to  just  l/8th  of  the  unit  cell  as  shown  in  Figure  9.4.  For  all 
the  results  discussed  in  this  section,  the  analysis  domain,  which  is  the  l/8th  unit  cell,  is 
part  of  the  bottom  ply  in  the  configuration.  Therefore,  the  bottom  surface  of  the  model  in 
Figure  9.4  is  traction-free  and  exposed  to  oxygen. 
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9.4.  Material  system 

The  material  system  used  for  all  the  analyses  discussed  in  this  section  is  Graphite/PMR- 
15  composite.  The  coupled  model  requires  both  the  oxidation  material  properties  as  well 
as  the  mechanical  properties  for  the  tow  and  matrix  in  the  composite.  Note  that  the 
configuration  is  assumed  to  be  at  a  temperature  of  288  C  throughout  the  entire 
simulation.  The  coupled  model  also  requires  the  degradation  schemes  for  the  matrix  and 
the  tow,  for  both  the  oxidation  as  well  as  mechanical  damage.  These  degradation 
schemes  are  described  in  the  next  section. 


Figure  9.4:  Analysis  domain  (l/8th  unit  cell)  with  transparent  matrix 


The  oxidation  material  properties  that  are  used  in  these  models  have  already  been 
described  in  the  previous  sections  that  discuss  the  oxidation  behavior  in  composites.  The 
oxidation  material  properties  for  the  neat  PMR-15  resin  were  obtained  from  ref  [25].  The 
oxidation  material  properties  for  the  tow  were  determined  using  the  homogenization 
strategies  described  in  Section  3.  Table  5.1  gives  the  oxidation  material  properties  for  the 
neat  PMR-15  resin  and  Table  7.1  gives  the  corresponding  properties  for  the  tow. 

Obtaining  the  mechanical  properties  for  the  Graphite/PMR-15  material  system  at  288  C 
is  not  easy  since  they  tend  to  change  over  time  and  it  may  not  be  appropriate  to  use 
property  data  from  different  sources  or  manufacturers  over  different  time  periods.  That 
being  the  case,  it  is  also  very  difficult  to  obtain  the  entire  set  of  required  mechanical 
properties  from  one  source  in  the  literature.  Moreover,  some  of  the  required  properties  at 
288  C  are  unavailable  due  to  the  lack  of  appropriate  experimental  techniques  to 
determine  them.  The  resin  and  the  tow  are  also  assumed  to  be  linear  elastic  materials 
although  elasto-plastic  behavior  of  the  polyimide  resin  would  be  expected  to  be  more 
prominent  at  288  C.  Overall,  the  set  of  mechanical  properties  for  the  material  system 
used  in  this  work  was  chosen  from  a  combination  of  different  sources  in  the  literature 
and  based  on  certain  assumptions  and  estimates  that  are  described  below. 
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The  mechanical  properties  of  the  neat  PMR-15  resin  were  chosen  based  on  experimental 
data  from  Pochiraju  and  Tandon  [30].  The  Young’s  modulus  of  the  neat  resin  was  found 
to  be  2.096  GPa  and  the  Poisson’s  ratio  is  assumed  to  be  0.30  in  Pochiraju  and  Tandon’s 
work  [30].  Based  on  the  assumption  that  the  neat  matrix  is  isotropic,  the  Young’s 
modulus  and  Poisson’s  ratio  can  be  used  to  calculate  the  shear  modulus.  Pochiraju  and 
Tandon  [30]  also  provide  the  normal  strength  at  room  temperature  and  288  C.  The  shear 
strength  of  the  neat  PMR-15  resin  is  calculated  by  scaling  the  strength  at  room 
temperature  based  on  the  change  in  normal  strength  from  room  temperature  to  288  C. 
Table  9.1  contains  the  elastic  moduli  for  neat  PMR-15  resin  that  were  used  in  this  work. 
The  strength  properties  that  were  discussed  in  this  paragraph  are  provided  under  Set  1  in 
Table  9.2.  The  properties  under  Set  2  and  the  need  for  an  additional  set  of  strength 
properties  are  discussed  in  the  next  paragraph. 

The  Graphite/PMR-15  tow  is  assumed  to  be  transversely  isotropic  and  therefore  its 
elastic  behavior  is  defined  by  five  independent  properties.  The  engineering  properties  for 
the  tow  were  harder  to  obtain  because  the  configuration  requires  properties  at  288  C.  The 
elastic  moduli  chosen  were  interpolated  from  work  performed  by  Odegard  and  Kumosa 
[50],  which  looked  at  the  effect  of  temperature  on  some  of  the  engineering  properties  of 
a  Graphite/PMR-15  unidirectional  laminate  (Vf=55.6%).  Of  the  five  independent 
properties  required,  En,  E22,  V12  and  G12  were  obtained  by  interpolating  from  the  data  in 
Ref  [50],  The  Poisson’s  ratio  in  the  transverse  plane,  V23,  at  288  C  was  assumed  to  be  the 
same  as  that  at  room  temperature.  The  elastic  moduli  for  the  tow  material  are 
summarized  in  Table  9.1.  All  the  strengths  properties  of  the  tow  at  288C  were  not 
available  in  the  literature.  It  is  relatively  difficult  to  determine  all  the  strength  properties 
for  the  tow.  These  properties,  especially  the  matrix-dominated  properties,  are  hard  to 
determine,  because  of  many  factors  like  the  material  interface  properties  that  influence 
the  strengths.  The  cr22  strength  and  <jn  strength  were  interpolated  from  Odegard  and 

Kumosa’ s  work[50].  Since,  the  tow  is  assumed  to  be  transversely  isotropic,  the  <r33 
strength  is  the  same  as  the  cr22  strength  and  the  cru  strength  is  the  same  as  the  cr13 
strength.  Due  to  lack  of  experimental  data  for  the  fiber-dominated  crn  strength,  the 

corresponding  strength  from  a  Graphite/epoxy  material  system  was  used.  However,  the 
<Tn  strength  is  only  consequential  only  during  fiber-breakage  which  occurs  during  final 

failure  of  the  composite.  Therefore,  this  assumption  was  not  considered  to  be  significant 
because  this  work  is  more  concerned  with  the  damage  initiation  and  progression  then  the 
final  failure  of  the  composite.  Also,  due  to  lack  of  experimental  data  for  the  <x23  strength, 

it  was  assumed  to  be  the  <x23  stress  corresponding  to  the  same  strain  level  at  which  the 
an  stress  mode  failed.  These  strength  properties  for  the  tow  are  summarized  under  Set  1 
in  Table  9.2.  Note  that  the  cr22  strength  of  the  tow  in  Set  1  is  considerably  lower  than  the 

normal  strength  of  the  neat  resin.  This  would  indicate  that  the  tows  would  fail  before  the 
neat  resin  pockets  in  the  composite.  It  is  common  for  a  composite  to  have  a  lower 
transverse  tensile  strength  than  the  tensile  strength  of  the  neat  resin  [51].  However,  since 
the  properties  for  this  material  system  were  compiled  from  different  sources  and 
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therefore  as  mentioned  before,  not  particularly  reliable,  another  set  of  assumed  strengths 
were  also  chosen  for  the  material  system.  In  this  new  set  of  properties,  the  normal  and 
shear  strength  of  the  neat  resin  were  scaled  down  based  on  typical  strength  ratios 
between  resin  and  tow  transverse  strengths  in  Graphite/Epoxy  material  systems.  This 
additional  set  of  strength  properties  for  the  material  system  used  in  this  work  is  defined 
as  Set  2  in  Table  9.2.  Having  two  sets  of  material  properties  would  also  give  another 
perspective  on  the  damage  initiation  and  progression  behavior  based  on  the  change  in 
engineering  properties. 


Table  9.1:  Elastic  properties  for  the  Graphite/PMR-15  material  system  [50,30] 


Resin 

Neat  PMR-15 

Tow 

Graphite/PMR- 1 5 

Ell 

2.096  GPa 

124.05  GPa 

E22=E33 

2.096  GPa 

6.2  GPa 

G12=G13 

0.806  GPa 

1.62  GPa 

G23 

0.806  GPa 

1.929  GPa 

vl2=vl3 

0.3 

0.485 

v23 

0.3 

0.607 

Table  9.2:  Strength  properties  for  the  Graphite/PMR-15  material  system  [50,30] 


Set  1 

Set  2 

Resin 

Neat  PMR-15 

Tow 

Graphite/PMR- 1 5 

Resin 

Neat  PMR-15 

Tow 

Graphite/PMR- 1 5 

Sll 

41 

2550 

12 

2550 

S22 

41 

18.91 

12 

19 

S33 

41 

18.91 

12 

19 

S12 

73.72 

36.83 

12 

37 

S23 

73.72 

43.85 

12 

44 

S13 

73.72 

36.83 

12 

37 

All  strengths  in  MPa 


9.5.  Constitutive  relations 

This  section  describes  the  different  constitutive  relations  that  are  required  to  implement 
this  coupled  oxidation-damage  progression  model.  This  includes  the  mechanical 
property  degradation  schemes  both  due  to  mechanical  loading  as  well  as  due  to  the  effect 
of  oxidation.  The  section  also  talks  about  how  the  two  degradation  schemes  are 
combined  and  used  in  the  coupled  model  to  obtain  the  overall  mechanical  properties  of 
the  material  based  on  the  oxidation  and  damage  state. 
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9.5.1.  Property  degradation  based  on  mechanical  damage 

When  a  failure  mode  is  triggered  during  the  damage  analysis,  the  engineering  properties 
are  degraded  to  account  for  the  change  in  mechanical  behavior.  This  operation  is 
conducted  based  on  a  property  degradation  scheme,  which  has  been  briefly  discussed  in 
Section  3.  In  this  work,  the  failure  criteria  and  the  property  degradation  scheme  is 
applied  on  each  integration  point  within  every  element  in  the  model.  The  property 
degradation  scheme  is  implemented  such  that  a  material  point  that  has  already  failed 
under  a  particular  mode  can  fail  under  another  mode.  In  such  a  case,  the  material 
properties  are  degraded  based  on  which  failure  mode  prescribes  the  larger  degradation. 

Different  property  degradation  schemes  have  been  proposed  in  the  literature  by  several 
researchers  such  as  Whitcomb  et  al.  [52],  Blackketter  et  al.  [53],  Tamma  et  al.  [54]  and 
Zako  et  al.  [55].  All  these  models  share  certain  similarities  and  differences.  They  are 
similar  in  the  sense  that  each  of  them  controls  the  amount  of  degradation  under  different 
failure  modes  for  both  the  tow  and  the  matrix.  Goyal  [49]  performed  a  comparison  of  the 
different  degradation  schemes  and  developed  a  common  framework  that  allowed 
implementation  of  a  wide  variety  of  degradation  schemes. 

For  all  the  damage  progression  models  in  this  work,  the  degradation  scheme  by 
Blackketter  [53]  was  used.  This  type  of  degradation  scheme  has  been  widely  used  by 
many  researchers  [52,56-58]  to  predict  initiation  and  progression  of  damage.  The 
engineering  properties  are  degraded  as  specified  in  Eq.(3.6).  The  degradation  scheme 
used  is  different  for  the  tow  and  the  matrix. 

In  the  degradation  scheme  for  the  tow  material,  the  values  of  the  degradation  parameters, 
cii  (i=l  to  6)  are  1,  5  or  100.  Note  that  the  value  of  the  parameters  in  a,  will  be  different 
under  different  damage  modes.  Table  9.3  gives  the  values  of  the  degradation  factors  for 
the  tow  material  under  this  scheme.  The  “1”  is  the  local  fiber  direction  of  the  tow  and 
“2”  and  “3”  are  the  local  transverse  directions  of  the  tow.  An  a,  value  of  1  indicates  that 
the  modulus  has  not  been  degraded.  An  a,  value  of  5  indicates  the  modulus  has  been 
degraded  to  20%  of  its  original  value  and  similarly  an  a,  value  of  100  indicates  the 
modulus  has  been  degraded  to  1%  of  its  original  value.  The  reason  that  some  of  the 
moduli  are  degraded  to  1%  of  the  original  rather  than  an  absolute  zero  is  to  avoid 
numerical  instabilities[53].  The  shear  moduli  were  not  reduced  to  less  than  20%  of  the 
original  value  under  mode  <x22  and  <r33  failure  because  it  is  assumed  that  some  shear 

stiffness  remains  due  to  frictional  resistance  still  present  on  the  failure  plane  [53]. 

For  the  matrix  material,  the  property  degradation  was  assumed  to  be  the  same  under  all 
the  failure  modes.  The  tensile  moduli  and  Poisson’s  ratios  of  the  matrix  are  reduced  to 
1%  of  its  original  value  whereas  the  shear  moduli  are  reduced  to  20%  of  its  original 
value.  The  matrix  is  therefore  assumed  to  become  anisotropic  after  failure.  Table  9.4 
gives  the  degradation  factors  for  the  matrix  material. 
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Table  9.3:  Degradation  parameters  (a,-)  for  engineering  elastic  properties  of  the  tow  [53] 


Mode  (7U 

Mode  cr22 

Mode  <t3 3 

Mode  <t12 

Mode  <t23 

Mode  <r13 

Ell 

100 

i 

i 

1 

i 

1 

E22 

100 

100 

i 

100 

100 

1 

E33 

100 

1 

100 

1 

100 

100 

G12 

100 

5 

1 

100 

100 

1 

G23 

100 

5 

5 

1 

100 

1 

G13 

100 

1 

5 

1 

100 

100 

vl2 

100 

1 

1 

1 

1 

1 

v23 

100 

100 

1 

100 

100 

1 

vl3 

100 

1 

1 

1 

1 

1 

Table  9.4:  Degradation  parameters  (a,)  for  engineering  elastic  properties  of  matrix  [53] 


All  Modes 

Ell 

100 

E22 

100 

E33 

100 

G12 

5 

G23 

5 

G13 

5 

vl2 

100 

v23 

100 

vl3 

100 

9.5.2.  Property  degradation  based  on  oxidation  damage 

The  effect  of  oxidation  on  the  mechanical  behavior  of  the  composites  is  considered  in 
the  coupled  analysis  models  used  in  this  work.  In  reality,  the  mechanical  behavior  is 
probably  more  tightly  coupled  with  the  oxidation  behavior  than  what  is  assumed  in  the 
current  model  because  the  mechanical  damage  can  affect  the  oxidation  behavior  by 
allowing  more  oxygen  to  penetrate  the  composite  material.  This  can  further  affect  the 
mechanical  behavior  because  more  oxidation  will  cause  more  damage  in  the  composite. 
These  complex  effects  are  not  considered  in  this  current  work.  In  this  work,  the 
oxidation  is  assumed  to  affect  the  mechanical  behavior,  but  not  the  converse. 

A  simple  constitutive  relation  or  property  degradation  scheme  was  developed  to  account 
for  the  effect  of  oxidation  on  the  mechanical  behavior  and  is  described  in  Section  3.  For 
a  general  orthotropic  material,  the  engineering  moduli  are  modified  according  to 
Eq.(3.6).  Unlike  the  property  degradation  scheme  for  mechanical  damage,  there  is  no 
failure  criteria  on  which  the  degradation  scheme  is  based. 

While  the  property  degradation  scheme  due  to  mechanical  damage  typically  reduces  the 
value  of  the  moduli,  the  same  is  not  necessarily  the  case  for  the  property  degradation 
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scheme  for  oxidation.  Experimental  work  has  shown  that  the  stiffness  of  the  fully 
oxidized  matrix  is  typically  larger  than  that  of  the  un-oxidized  material  [25].  There  is  not 
enough  data  in  the  literature  in  order  to  determine  all  the  degradation  parameters,  bj.  In 
order  to  implement  the  degradation  scheme  for  this  work,  the  best  available  data  in  the 
literature  was  used  where  appropriate  and  estimates  based  on  certain  assumptions  were 
used  to  the  remaining  parameters.  The  values  of  bj  chosen  for  the  matrix  and  tow 
materials  in  this  work  are  given  in  Table  9.5.  Experiments  showed  that  the  elastic 
modulus  of  the  neat  PMR-15  resin  increased  by  about  20%  when  fully  oxidized  [25], 
The  same  amount  of  increase  is  assumed  to  apply  for  the  shear  moduli.  The  Poisson’s 
ratio  is  assumed  to  remain  constant  based  on  the  assumption  that  the  matrix  remains 
isotropic  after  oxidation.  The  same  challenges  exist  for  obtaining  accurate 
characterization  data  for  tows  or  unidirectional  laminates.  The  fiber  is  assumed  to  be 
impermeable  and  unaffected  by  the  oxidation.  Simple  micromechanics  analyses  showed 
that  effective  tow  properties  were  changed  by  a  very  small  amount  when  the  matrix 
moduli  were  increased  by  20%.  Since  the  change  was  insignificant,  the  degradation 
properties  (b,)  for  the  tow  were  assumed  to  be  zero,  meaning  that  the  elastic  properties  of 
the  tow  were  assumed  to  remain  constant  after  oxidation.  Therefore,  an  undamaged 
material  point  in  the  tow  was  assumed  to  remain  transversely  isotropic  after  the  material 
was  oxidized. 

As  mentioned  in  Section  3,  the  property  degradation  scheme  based  on  oxidation  also 
degrades  the  strength  properties  of  the  materials  in  the  composite  as  defined  by 
Eq.(3.103).  There  is  no  data  in  the  literature  that  can  be  used  to  determine  the  strength 
degradation  parameters,  dj.  Due  to  this  limitation,  for  all  the  models  analyzed  in  this 
work,  strength  properties  are  assumed  for  the  fully  oxidized  matrix  and  tow.  Table  9.6 
gives  the  values  of  the  strength  degradation  parameters  chosen  for  the  matrix  and  tow. 
The  strengths  for  all  stress  components  in  the  matrix  are  assumed  to  drop  to  half  its 
value.  In  the  case  of  the  tow  material,  the  crn  strength,  which  is  the  strength  in  the  fiber 

direction,  is  assumed  to  drop  to  95%  of  the  original  value  whereas  all  the  other  strengths 
drop  50%. 
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Table  9.5:  Parameters  ( bi )  for  degrading  engineering  elastic  properties  of  matrix  and  tow 


Engineering  property 
affected 

bi 

Matrix 

Tow 

1 

Ell 

+0.2 

0.0 

2 

E22 

+0.2 

0.0 

3 

E33 

+0.2 

0.0 

4 

G12 

+0.2 

0.0 

5 

G23 

+0.2 

0.0 

6 

G13 

+0.2 

0.0 

7 

vl2 

0 

0.0 

8 

v23 

0 

0.0 

9 

vl3 

0 

0.0 

Table  9.6:  Parameters  (di)  for  degrading  strength  properties  of  matrix  and  tow 


Strength  property 

di 

affected 

Matrix 

Tow 

1 

Sll 

-0.50 

-0.05 

2 

S22 

-0.50 

-0.50 

3 

S33 

-0.50 

-0.50 

4 

S12 

-0.50 

-0.50 

5 

S23 

-0.50 

-0.50 

6 

S13 

-0.50 

-0.50 

The  overall  mechanical  moduli  of  the  material  are  obtained  by  combining  the  two 
degradation  schemes,  both  based  on  mechanical  damage  as  well  as  oxidation,  as 
described  in  Section  3. The  expressions  for  the  overall  properties  at  a  material  point  are 
given  by  the  Eq.(3.104).  Note  that  although  the  degradation  scheme  chosen  in  this  work 
assumes  that  the  matrix  remains  isotropic  after  oxidation,  the  overall  mechanical 
properties  obtained  after  accounting  for  mechanical  damage  need  not  necessarily 
represent  an  isotropic  material.  The  parameters,  a„  have  a  value  of  1  if  the  material  is  not 
damaged  and  therefore  in  such  a  case,  the  matrix  would  remain  isotropic.  On  the  other 
hand,  if  the  matrix  is  damaged  under  any  mechanical  failure  mode,  the  matrix  becomes 
anisotropic.  Similarly,  the  tow  need  not  remain  transversely  isotropic  after  the 
mechanical  properties  have  been  modified  using  Eq.(3.104). 

9.6.  Results  and  discussion 

The  coupled  analysis  model  was  used  to  simulate  damage  initiation  and  progression  in 
the  configuration  described  in  Section  9.3.  The  basic  configuration  described  in  Section 
9.3  is  a  two-ply  laminate  at  288  C  with  the  top  and  bottom  surfaces  exposed  to  oxygen. 
The  laminate  is  assumed  to  be  infinite  in  the  in-plane  directions  and  has  a  uniaxial  load 
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in  the  x-direction.  As  described  in  Section  9.4,  two  sets  of  material  properties  were 
chosen  to  model  the  Graphite/PMR-15  material  system.  This  section  will  discuss  the 
results  from  the  analyses  performed  using  both  sets  of  properties.  A  parametric  study 
was  also  performed  where  the  number  of  plies  in  the  laminate  was  increased.  The 
parametric  study  looked  at  two-,  four-  and  six-ply  laminates  for  both  the  sets  of  material 
properties.  The  results  of  this  parametric  study  will  be  described  in  this  section  as  well. 


9.6.1.  Two-ply  laminate 

The  damage  progression  behavior  of  the  laminate  under  mechanical  load  alone  (i.e.  no 
oxidation)  is  first  discussed.  The  laminate  is  assumed  to  be  quasi-statically  loaded 
uniaxially  while  maintained  at  a  temperature  of  288  C.  Since  two  sets  of  material 
properties  were  chosen  to  define  the  Graphite/PMR-15  material  system,  the  damage 
analyses  were  performed  on  two  models,  one  for  each  material  property  set.  Note  that 
the  two  sets  of  material  properties  have  the  same  elastic  moduli.  The  difference  between 
the  two  sets  of  material  properties  is  in  the  strengths  properties  as  shown  in  Table  9.2. 
Figure  9.5  shows  a  plot  of  the  volume  averaged  <yxx  versus  the  volume  averaged  sxx  for 

both  the  models.  As  expected,  the  stress-strain  behaviors  are  different  for  the  two 
models.  Figure  9.5  shows  that  the  initial  damage  in  the  model  using  Set  1  properties 
causes  a  significant  drop  in  load  (indicated  by  A)  compared  to  the  initial  damage  in  the 
model  using  Set  2  properties  (indicated  by  B).  This  difference  in  behavior  can  be 
explained  by  looking  at  where  the  initial  damage  occurs.  In  the  case  of  Set  1,  in  which 
the  transverse  tow  strengths  are  much  lower  than  the  matrix  strengths,  damage  initiated 
in  the  fill  tow  under  compressive  cr33  damage  mode.  The  observation  that  parts  of  the 

fill  tow  closer  to  the  laminate  mid-plane  are  under  compression  can  be  explained  by 
considering  that  warp  tows  are  being  stretched  because  of  the  load  and  therefore  pushing 
on  the  fill  tows  in  between.  When  the  material  properties  of  the  damaged  area  in  the  fill 
tow  are  degraded,  the  amount  of  load  carried  by  the  tow  reduces.  In  the  case  of  Set  2,  in 
which  the  resin  has  the  lowest  strengths,  the  damage  initiates  in  the  matrix  pockets  under 
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Figure  9.5:  Volume  averaged  stress-volume  averaged  strain  curves  for  the  two- 

ply  laminate  without  oxidation 


tensile  cr33  damage  mode  but  since  the  matrix  doesn’t  carry  as  much  load  as  the  tows, 

the  load  drop  is  not  as  significant  as  that  seen  in  the  model  with  Set  1  properties.  This 
explanation  can  be  confirmed  by  looking  at  the  damage  evolution  in  the  two  models. 

Figure  9.6  shows  the  evolution  of  damage  in  the  different  constituents  of  the  model 
using  Set  1  material  properties.  It  shows  the  location  of  the  initial  damage  in  the  fill  tows 
at  a  volume  averaged  sxx  strain  level  of  0.0935%  strain.  The  initial  damage  occurs  under 

compressive  cr33  damage  mode.  It  can  be  seen  that  the  matrix  is  the  last  constituent  in 

the  composite  to  have  significant  failure.  Looking  at  the  column  for  0.6%  strain  in 
Figure  9.6  shows  that  the  there  is  significant  transverse  damage  in  the  fill  and  warp  tows 
but  there  is  hardly  any  damage  in  the  matrix.  This  behavior  was  also  expected  based  on 
the  fact  that  the  Set  1  material  properties  have  the  transverse  tow  strengths  much  lower 
than  that  of  the  matrix. 

In  comparison,  Figure  9.7  shows  the  evolution  of  damage  in  the  model  using  Set  2 
material  properties.  In  this  case,  it  shows  that  the  damage  initiates  in  the  matrix  under 
tensile  cr33  mode  near  the  mid-plane  of  the  laminate  at  a  volume  averaged  sxx  strain 

level  of  0.0473%.  Note  that  the  damage  initiates  at  a  much  lower  strain  level  when  using 
Set  2  material  properties  versus  those  in  Set  1.  Although  the  initial  damage  in  the  Set  2 
model  is  in  the  matrix,  the  first  significant  drop  in  load  is  at  a  strain  level  of  0.128% 
(indicated  by  C  in  Figure  9.5)  and  it  is  caused  by  damage  in  the  fill  tow  under 
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compressive  cr33  failure  mode.  In  comparison,  the  first  significant  drop  in  the  Set  1 

model  occurs  at  0.0935%  strain  (indicated  by  A  in  Figure  9.5)  under  the  same  type  of 
failure  mode.  The  reason  why  the  damage  in  the  fill  tow  occurs  at  a  higher  strain  level  in 
the  Set  2  model  could  be  explained  as  follows.  When  the  Set  2  laminate  is  at  0.0935% 
strain,  there  is  already  some  damage  in  the  matrix  pockets.  This  would  make  the  matrix 
pockets  more  compliant  and  thereby  effectively  reducing  the  constraints  on  the  fill  tow. 
The  fill  tows  would  be  allowed  to  deform  more  freely  than  before  matrix  damage 
occurred  and  therefore  relieving  the  cr33  stresses  in  the  fill  tow.  Thus,  a  larger  load 

would  be  required  to  raise  the  cr33  stress  in  the  fill  tow  enough  to  cause  damage. 

In  comparison  to  the  evolution  of  damage  in  the  Set  1  model,  the  matrix  has  much  more 
damage  at  0.6%  strain.  Also,  there  are  slight  differences  between  the  damaged  locations 
in  the  tows.  This  is  probably  because  of  the  manner  in  which  the  load  is  transferred  when 
different  locations  in  the  laminate  start  to  fail. 

The  results  discussed  up  to  this  point  considered  only  the  effect  of  damage  due  to 
mechanical  loading.  Now  the  results  from  the  coupled  models  are  discussed.  The 
coupled  models  simulate  the  mechanical  behavior  when  the  laminate  is  under  a  fixed 
mechanical  loading  and  is  then  exposed  to  oxygen  from  the  top  and  bottom  surfaces  for 
200  hours.  These  simulations  are  performed  at  different  fixed  mechanical  loads.  Similar 
simulations  are  performed  on  models  with  each  set  of  material  properties  to  determine 
the  effect  of  the  properties  on  the  behavior. 

The  behavior  of  the  model  using  Set  1  material  properties  is  discussed  first.  As 
illustrated  in  Figure  9.6,  damage  due  to  a  mechanical-only  load  initiates  at  a  strain  level 
of  0.0935%.  A  coupled  model  simulation  was  performed  at  a  strain  level  of  0.09%  to  see 
if  the  oxidization  would  initiate  any  damage.  It  was  seen  that  there  was  no  effect  of 
oxidation  on  the  damage  behavior  throughout  the  200  hours.  This  is  because,  as  shown 
in  Figure  9.6,  all  the  initial  damage  is  located  in  the  top  half  of  the  fill  tow  in  the  model, 
which  implies  that  the  stress  failure  index  is  highest  in  that  region  of  the  fill  tow.  This 
region  in  the  model  corresponds  to  the  interior  of  the  laminate  because  the  analysis 
domain  represents  the  lower  half  of  the  laminate.  After  200  hours  of  oxidation,  the 
oxidation  front  has  not  reached  the  interior  of  the  laminate  far  enough  to  affect  the 
engineering  properties  of  the  tow  to  cause  damage.  As  defined  in  Table  9.5,  the  change 
in  engineering  moduli  is  not  significant  enough  to  affect  the  stresses.  The  changes  in  the 
strength  properties  are  significant,  but  the  regions  with  the  stress  concentrations  are 
either  not  oxidized,  or  not  oxidized  enough  to  cause  damage  in  the  fill  tows. 
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Strain  Level 


0.0935  %  |  0.2  % 


0.4% 


0.6% 


Fill 

Tow 


Warp 

Tow 


Failure  mode  33 


Figure  9.6:  Evolution  of  damage  in  the  two-ply  laminate  configuration  without  oxidation  using  Set  1  material  properties 


Failure  mode  33 


Warp 

Tow 


Figure  9.7:  Evolution  of  damage  in  the  two-ply  laminate  configuration  without  oxidation  using  Set  2  material  properties 
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The  simulations  were  also  performed  at  strain  levels  of  0.2%  and  0.4%.  Since  the 
configuration  is  assumed  to  be  already  loaded  to  a  constant  strain  level,  the  configuration 
should  also  be  assumed  to  have  the  damage  that  would  have  ordinarily  occurred  without 
the  influence  of  oxidation.  This  initial  damage  state  for  a  particular  load  level  is  assumed 
based  on  the  corresponding  damage  state  for  the  model  from  the  standard  damage 
progression  analysis.  In  the  model  with  Set  1  material  properties  at  a  strain  level  of 
0.2%,  the  damage  at  the  beginning  of  the  oxidation  simulation  is  almost  entirely  in  the 
fill  tow  as  shown  in  Figure  9.6.  There  is  no  damage  at  all  in  the  matrix.  The  only  other 
damage  in  the  configuration  is  one  integration  point  in  each  warp  tow  that  has  failed  in 
the  cr33  damage  mode  as  shown  in  Figure  9.6.  Figure  9.8  shows  the  evolution  of  damage 

as  the  oxidation  progresses  for  200  hours.  After  one  hour  of  oxidation,  there  is  new 
damage  under  cr22  and  a33  failure  modes  in  the  bottom  region  of  the  fill  tow  where  the 

oxygen  is  slowly  making  its  way  into  the  interior  of  the  laminate.  The  simulation  also 
shows  some  slight  damage  in  the  matrix  pocket  closer  to  the  exposed  surface  of  the 
laminate.  There  is  also  some  damage  in  the  <x22  failure  mode  in  the  lower  half  of  the 

warp  tow,  which  can  be  explained  due  to  the  oxidation  front  creeping  into  the  interior  of 
the  laminate.  The  more  interesting  behavior  is  that  regions  of  the  top  half  of  the  warp 
tow  fails  in  the  <j33  damage  mode.  This  is  interesting  because  the  damage  is  seen  after 

only  one  hour  of  oxidation,  at  which  time  the  oxidation  front  has  not  reached  even  close 
to  the  top  half  of  the  model.  This  can  be  explained  by  the  redistribution  of  the  load  in  the 
configuration  after  material  damage.  As  mentioned  earlier,  even  before  oxidation  began, 
there  was  damage  in  the  fill  tow.  Figure  9.8  shows  that  after  only  one  hour  of  oxidation, 
there  is  significant  damage  in  the  fill  tow,  which  renders  most  of  the  fill  tow  incapable  of 
carrying  load.  This  increases  the  load  on  the  warp  tow.  The  effect  of  the  external  load  on 
the  laminate  is  to  straighten  the  undulating  warp  tows,  which  causes  a  tensile  cr33  in  the 
top  half  of  the  tow.  When  the  load  on  the  warp  tow  increases,  it  also  increases  the  <j33 

stress  in  the  top  half  of  the  tow  making  it  exceed  the  strength.  This  behavior  shows  that 
the  influence  of  oxidation  on  the  mechanical  behavior  is  not  always  localized  and  in 
some  cases,  its  effect  can  be  seen  in  the  interior  of  the  laminate  where  the  material  has 
not  been  oxidized. 

Figure  9.9  shows  the  initial  damage  state  in  the  two-ply  laminate  at  0.4%  strain  before 
oxidation  begins.  It  shows  that  there  is  very  little  damage  in  the  matrix  pockets.  The  fill 
tow  on  the  other  hand  has  considerable  damage  in  the  cr22  and  <j33  failure  modes.  Figure 

9.9  also  shows  that  warp  tow  has  some  damage  in  the  top  half  under  mainly  the  <r33 
failure  mode.  As  the  oxidation  progresses,  some  build-up  of  uu  damage  is  seen  in  the 

lower  matrix  pocket  as  shown  in  Figure  9.9.  There  is  little  new  damage  in  the  fill  tow 
since  most  of  the  tow  was  already  damaged  before  the  oxidation  began.  The  warp  tow 
sees  considerable  new  damage  under  the  cr22  failure  mode  in  the  bottom  half  of  the  tow 

as  oxidation  progresses.  This  can  be  explained  by  the  fact  that  the  fill  tow  is  mostly 
damaged  and  much  of  the  load  is  now  carried  by  the  warp  tow.  Therefore,  the  warp  tow 
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would  experience  higher  stresses.  In  addition  to  the  higher  stresses,  the  oxidation  causes 
the  strengths  to  drop  by  50%  thereby  increasing  the  potential  of  failure. 

Figure  9.10  shows  the  plot  of  the  volume  averaged  <rxx  with  the  volume  averaged  sxx  for 
the  Set  1  model  indicating  the  drop  in  the  volume  average  axx  at  the  end  of  200  hours  of 

oxidation  for  the  two  simulations  discussed  earlier.  It  shows  that  for  the  0.2%  strain  level 
simulation,  the  volume  averaged  crxx  stress  reduced  from  the  point  labeled  A  to  A’ 

indicating  a  drop  of  15%  at  the  end  of  200  hours  of  oxidation  whereas  in  the  case  of  the 
0.4%  strain  level,  the  corresponding  stress  dropped  over  22%  indicated  by  the  line  B-B’. 
Figure  9.11  shows  the  volume  averaged  <jxx  for  all  three  simulations  normalized  with  the 

initial  volume  averaged  stress  as  the  oxidation  progresses  over  200  hours.  As  mentioned 
earlier,  at  the  0.09%  strain  level,  there  was  no  new  damage  due  to  oxidation  and 
therefore  there  was  no  drop  in  the  volume  averaged  stress.  Instead,  there  was  a  slight 
increase  in  the  volume  averaged  stress  due  to  the  fact  that  the  stiffness  in  the  matrix 
increases  when  oxidized  but  the  increase  is  so  small  that  it  is  not  noticeable  in  Figure 
9.11.  In  the  case  of  the  0.2%  strain  level,  a  significant  part  of  the  stress  drop  occurs  in 
the  beginning  of  the  oxidation  process  within  the  first  two  hours.  This  indicates  that  the 
damage  that  occurred  in  the  remaining  198  hours  was  not  significant  enough  to  reduce 
the  load  in  the  laminate.  In  the  case  of  the  0.4%  strain  level,  a  major  portion  of  the  stress 
drop  occurs  at  a  single  time  step  at  53.33  hours  when  the  stress  drops  to  78.3%  of  the 
initial  stress.  The  damage  that  occurs  before  and  after  that  point  accounts  for  just  0.7% 
of  the  total  drop  in  load. 

The  results  from  the  simulations  of  the  laminate  using  the  Set  2  material  properties  are 
discussed  next.  Similar  to  the  simulations  on  the  laminates  with  Set  1  material 
properties,  three  simulations  were  performed  with  strain  levels  of  0.1%,  0.2%  and  0.4%. 
In  the  simulation  with  0.1%  strain,  the  initial  damage,  as  shown  in  Figure  9.12,  is 
confined  to  mostly  the  inter-laminar  matrix  pocket.  There  is  also  slight  damage  under 
<r33  failure  mode  in  the  top  part  of  the  warp  tow.  The  coupled  analysis  shows  that  there 

is  no  new  damage  caused  due  to  the  effect  of  oxidation.  This  is  similar  to  the  model  with 
Set  1  material  properties  and  0.09%  strain,  where  the  stress  state  in  the  oxidized 
materials  is  not  significant  enough  to  cause  new  damage.  With  0.2%  strain,  the  initial 
damage  is,  as  expected,  more  widespread  than  that  in  the  case  with  0.1%  strain.  As 
shown  in  Figure  9.13,  the  damage  in  the  inter-laminar  matrix  pocket  has  increased  in 
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Figure  9.9:  Evolution  of  damage  due  to  oxidation  in  the  two-ply  laminate  at  0.4%  strain  using  Set  1  material  properties 
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Figure  9.10:  Volume  averaged  stress-volume  averaged  strain  for  the  Set  1 
material  two-ply  laminate  showing  drop  in  stress  after  200  hours  of  oxidation 


Figure  9.1 1:  Variation  in  volume  averaged  stress  due  to  oxidation  for  the 
Set  1  material  two-ply  laminate  at  different  strain  levels 
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addition  to  new  <yn  failure  in  the  bottom  matrix  pocket.  The  fill  tow  has  <x22  and  <r33 

failure  in  regions  from  the  bottom  to  the  top,  with  more  damage  in  the  latter.  Figure  9.13 
also  shows  that  the  damage  in  the  warp  tow  increased.  The  simulation  predicted  that  the 
damage  in  the  lower  matrix  pocket  grows  mostly  under  <jn  and  cr22  failure  modes.  The 

simulation  also  predicts,  as  shown  in  Figure  9.13,  that  there  is  new  damage  in  the  bottom 
part  of  the  fill  tow  under  the  <j22  failure  mode.  This  can  be  explained  as  a  direct  effect 

of  the  oxidation  of  the  tows  which  reduces  the  strength  by  as  much  as  50%.  The  warp 
tow  also  has  new  damage  growth  as  an  effect  of  the  oxidation.  As  shown  in  Figure  9.13, 
the  warp  tow  starts  to  see  damage  in  the  lower  part  of  the  tow  under  cr22  failure  mode  as 

the  oxidation  progresses.  The  warp  tow  also  starts  to  have  <x33  damage  at  the  location 

indicated  by  A  as  the  oxidation  simulation  nears  the  end  of  200  hours.  A  similar 
simulation  was  performed  for  a  constant  strain  level  of  0.4%.  In  this  case,  the  initial 
damage  state  is  more  extensive  compared  to  the  simulation  with  0.2%  strain.  The  matrix, 
fill  tow  and  warp  tow  have  considerable  damage  as  shown  in  Figure  9.14.  As  expected, 
the  coupled  analysis  predicted  growth  in  the  damage  in  the  lower  matrix  pocket  due  to 
oxidation.  Most  of  this  damage  occurs  under  <ju  failure  mode  along  with  <x22  and  <r33 

failure  modes.  In  the  case  of  the  fill  tows,  the  oxidation  causes  additional  damage  in  the 
lower  part  of  the  tow  under  cr22  failure  mode,  as  shown  in  Figure  9.14.  Additional 

damage  is  also  seen  in  the  warp  tow  as  an  effect  of  the  oxidation.  The  new  damage  in  the 
warp  tow  occurs  under  cr22 ,  cr33  and  cr13  failure  modes. 

Similar  to  Figure  9.10,  Figure  9.15  shows  the  plot  of  the  volume  averaged  crxx  versus 
volume  averaged  sxx  for  the  Set  2  model  indicating  the  drop  in  the  volume  average  axx 

at  the  end  of  200  hours  of  oxidation  for  the  two  simulations  discussed  earlier.  Line  A- A’ 
shows  that  for  the  0.2%  strain  level  simulation,  the  volume  averaged  axx  stress  dropped 

13%  at  the  end  of  200  hours  of  oxidation  whereas  in  the  case  of  the  0.4%  strain  level,  the 
corresponding  stress  dropped  over  13.3%  (indicated  by  line  B-B’).  Figure  9.16  shows 
the  volume  averaged  axx  for  all  three  simulations  normalized  with  the  initial  volume 

averaged  stress  as  the  oxidation  progresses  over  200  hours.  Just  as  the  model  with  Set  1 
material  properties  at  the  0.09%  strain  level,  there  was  no  new  damage  due  to  oxidation 
and  the  volume  averaged  stress  actually  increases  slightly,  although  it  is  not  noticeable 
in  Figure  9.16.  In  the  case  of  the  other  two  strain  levels,  the  drop  in  volume  averaged 
stress  is  more  gradual  than  the  behavior  seen  in  the  corresponding  models  with  Set  1 
material  properties.  Although  there  are  some  sudden  drops  in  the  volume  averaged  stress 
as  seen  in  Figure  9.16,  they  are  not  as  significant  as  the  drops  seen  in  Figure  9.11.  This 
behavior  is  attributed  to  growth  in  the  matrix  damage  observed  in  the  Set  2  laminates 
that  is  not  seen  in  the  Set  1  laminates. 


Figure  9.12:  Evolution  of  damage  due  to  oxidation  in  the  two-ply  laminate  at  0.1%  strain  using  Set  2  material  properties 
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Figure  9.13:  Evolution  of  damage  due  to  oxidation  in  the  two-ply  laminate  at  0.2%  strain  using  Set  2  material  properties 
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Figure  9.14:  Evolution  of  damage  due  to  oxidation  in  the  two-ply  laminate  at  0.4%  strain  using  Set  2  material  properties 
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Figure  9.15:  Volume  averaged  stress-volume  averaged  strain  for  the  Set  2 
material  two-  ply  laminate  showing  drop  in  stress  after  200  hours  of 

oxidation 


Figure  9.16:  Variation  in  volume  averaged  stress  due  to  oxidation  for  the  Set  2 
material  two-ply  laminate  at  different  strain  levels 
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9.6.2.  Effect  of  number  of  plies 

A  parametric  study  was  conducted  to  see  if  the  number  of  plies  in  the  laminate  had  any 
significant  effect  on  the  mechanical  behavior  under  oxidation.  In  addition  to  the  laminates  that 
were  discussed  in  the  previous  section,  4-ply  and  6-ply  laminates  were  analyzed  for  both  Set  1 
and  Set  2  material  properties.  First,  a  damage  progression  analysis  (i.e.  no  oxidation)  was 
performed  on  the  laminates  to  determine  the  mechanical  behavior  as  the  load  on  the  laminate 
was  increased.  Figure  9.17  gives  a  plot  of  the  volume  averaged  crxx  with  the  volume  averaged 

sxx  for  the  different  laminates  using  Set  1  material  properties.  It  shows  that  the  overall  behavior 

is  not  significantly  different,  which  is  not  surprising.  Figure  9.18  shows  the  same  plot  for  the 
laminates  with  Set  2  material  properties.  Again,  the  number  of  plies  does  not  seem  to  have  an 
effect  on  the  overall  behavior.  Looking  at  the  evolution  of  damage  in  the  laminates  revealed 
generally  the  same  trends  as  seen  in  the  2-ply  laminates.  In  the  case  of  the  Set  1  laminates,  the 
fill  tows  had  initial  damage  and  continued  to  accumulate  much  more  damage  than  the  warp  tows 
followed  by  the  matrix,  which  had  very  little  failure.  In  the  case  of  Set  2  laminates,  as  seen  in  the 
corresponding  2-ply  laminate,  the  damage  initiates  in  the  matrix  followed  by  the  fill  tow  failing 
considerably  while  the  warp  tow  has  less  damage  in  comparison. 

The  coupled  simulations  were  performed  on  these  laminates  as  was  done  for  the  2-ply  laminates 
discussed  in  the  previous  section.  The  laminates  were  analyzed  at  different  strain  levels  and 
overall  they  showed  the  same  trends  as  seen  in  the  2-ply  laminates.  If  the  strain  levels  are  too 
low,  for  example  at  0.1%,  the  oxidation  was  not  found  to  have  any  significant  effect  of  the 
mechanical  behavior.  The  results  from  the  0.2%  and  0.4%  strain  level  simulations  will  be 
discussed  here.  Since  the  general  trends  are  the  same  as  compared  to  the  2-ply  laminates,  the 
evolution  of  damage  in  each  laminate  will  not  be  discussed  here.  Instead,  comparisons  of  the 
overall  behavior  will  be  discussed.  Comparing  the  results  from  the  2-ply,  4-ply  and  6-ply 
laminates  is  not  easy  since  they  do  not  follow  the  same  load  path  as  shown  in  Figure  9.18.  It 
would  definitely  not  make  sense  to  make  comparisons  at  same  strain  level  using  the  predicted 
volume  average  stress  values  because  of  the  same  reason.  It  makes  more  sense  to  look  at  the 
percentage  drop  in  the  volume  average  stress. 
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Figure  9.17:  Volume  averaged  stress-volume  averaged  strain  curves  for  the 
laminate  with  Set  1  material  properties 


Figure  9.18:  Volume  averaged  stress-volume  averaged  strain  curves  for 
the  laminate  with  Set  2  material  properties 
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Figure  9.19:  Variation  in  volume  averaged  stress  due  to  oxidation  for  the 
various  Set  1  material  laminates  at  0.2%  strain  level 


Figure  9.19  shows  the  variation  in  normalized  volume  average  stress  at  0.2%  strain  for  the  2-ply, 
4-ply  and  6-ply  laminates  using  Set  1  materials.  In  the  2-ply  laminate,  under  a  0.2%  strain  level, 
the  volume  average  stress  drops  15%  whereas  in  the  4-ply  model,  it  drops  only  3.4%,  which 
indicates  that  the  damage  in  the  4-ply  laminate  was  not  significant.  This  is  intuitive  since  a 
smaller  fraction  of  the  laminate  is  getting  oxidized  when  the  number  of  plies  increases  from  2  to 
4.  However  when  the  number  of  plies  is  increased  from  4  to  6,  the  stress  drop  increases  slightly 
from  3.4%  to  5.2%.  This  particular  trend  could  not  be  explained  but  as  discussed  later  in  this 
section,  this  counter-intuitive  behavior  was  not  observed  for  the  other  material  set  or  for  other 
strain  levels.  Figure  9.19  also  shows  that  all  the  Set  1  material  laminates  experience  the 
significant  drop  in  the  volume  average  stress  within  3  hours  of  oxidation.  Figure  9.20  shows  the 
variation  of  the  normalized  stress  for  the  0.4%  strain  level.  It  shows  the  percentage  drop  in  the 
volume  average  stress  at  the  end  of  200  hours  steadily  reducing  as  the  number  of  plies  in  the 
laminate  increase.  Comparison  of  Figures  9.19  and  9.20  shows  that  when  the  strain  level  was 
increased,  the  decrease  in  percentage  load  drop  was  more  gradual  with  the  number  of  plies.  In 
the  case  of  0.4%  strain,  the  drop  is  22.4%  for  a  2-ply  laminate,  14.8%  for  a  4-ply  laminate  and 
9.5%  for  a  6-ply  laminate.  On  the  other  hand,  in  the  0.2%  strain  level,  as  seen  in  Figure  9.19,  the 
percentage  drop  reduces  from  15%  to  less  than  6%  as  the  number  of  plies  is  increased  to  4  and  6. 
Note  that  this  trend  is  specific  to  the  laminates  with  Set  1  material  properties  and  cannot  be 
generalized.  Similar  simulations  were  performed  on  the  corresponding  laminates  with  Set  2 
material  properties  and  Figures  9.21  and  9.22  shows  the  variation  in  volume  average  stress  for 
0.2%  and  0.4%  strain  loading  respectively.  Again,  a  similar  trend  is  seen  where  there  is  a 
significant  reduction  in  the  percentage  drop  in  volume  average  stress  at  the  end  of  200  hours 
(from  13%  to  4.1%)  for  a  0.2%  loading  when  the  number  of  plies  is  increased  from  2  to  4.  When 
the  number  of  plies  is  increased  to  6,  the  drop  is  only  3.3%,  which  is  a  further  reduction  in  the 
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drop  compared  to  the  corresponding  case  in  the  Set  1  material  laminates.  Also,  the  variation  in 
the  normalized  stress  with  respect  to  time  is  much  more  gradual  in  the  Set  2  laminates  as 
compared  to  the  Set  1  laminates.  When  the  strain  level  is  increased  to  0.4%,  again  similar  trends 
are  seen  where  the  reduction  in  the  percentage  drop  is  not  as  drastic  as  that  for  0.2%  strain. 
Figure  9.22  shows  that  the  percentage  volume  average  stress  drop  at  the  end  of  200  hours  of 
oxidation  reduces  from  13.3%  to  7%  when  the  number  of  plies  go  up  from  2  to  4  and  the  drop 
further  reduces  to  5.2%  when  the  number  of  plies  is  increased  to  6.  This  behavior  can  be 
explained  based  on  the  fact  that  the  oxidation  process  oxidizes  the  same  amount  of  material  in  all 
these  laminates.  In  the  coupled  simulations  described  in  this  work,  the  oxidation  analysis  does 
not  depend  on  the  stress  or  damage  state  in  the  laminate.  Therefore,  regardless  of  the  number  of 
plies,  the  oxidation  layer  thickness  varies  in  the  same  manner  in  all  the  laminates.  This  also  has 
to  do  with  the  fact  that  at  the  end  of  the  200  hour  simulation,  the  maximum  predicted  oxidation 
layer  thickness  is  less  than  the  thickness  of  a  single  ply.  Increasing  the  number  of  plies  in  the 
laminate  effectively  increases  the  amount  of  material  that  can  carry  load,  but  the  oxidation 
process  only  affects  the  same  amount  of  material  regardless  the  number  of  plies.  Therefore,  it 
would  be  expected  that  the  percentage  drop  in  volume  average  stress,  or  load  drop,  would 
decrease  as  the  number  of  plies  increased. 
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Figure  9.20:  Variation  in  volume  averaged  stress  due  to  oxidation  for  the  various 
Set  1  material  laminates  at  0.4%  strain  level 
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Figure  9.21:  Variation  in  volume  averaged  stress  due  to  oxidation  for  the 
various  Set  2  material  laminates  at  0.2%  strain  level 


Figure  9.22:  Variation  in  volume  averaged  stress  due  to  oxidation  for  the 
various  Set  2  material  laminates  at  0.4%  strain  level 
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8.10.  Summary 

The  coupled  analysis  model  described  in  Section  3  is  used  to  predict  the  mechanical  behavior  of 
woven  composite  laminates  that  are  under  mechanical  load  as  well  as  exposed  to  oxygen.  The 
configuration  that  is  analyzed  and  the  complete  parameters  for  the  material  system  and  the 
constitutive  relations  are  described  in  this  section.  The  current  implementation  of  the  coupled 
analysis  model  makes  a  number  of  assumptions  when  simulating  the  behavior  of  the  laminate. 
The  effects  of  thermal  expansion  and  the  shrinkage  of  the  matrix  due  to  oxidation  are  ignored. 
These  are  effects  that  need  to  be  considered  in  future  implementations  of  the  coupled  analysis 
model  in  order  to  represent  more  accurately  the  behavior  of  the  underlying  mechanisms.  The 
effect  of  the  stress  and  damage  state  on  the  oxidation  behavior  also  needs  to  be  considered  in 
future  enhancements  of  the  coupled  model.  However,  the  analyses  described  in  this  work  provide 
a  framework  for  the  implementation  of  a  more  robust  tool  to  predict  the  behavior  of  laminates 
under  oxidation.  Due  to  lack  of  a  full  set  of  reliable  material  properties,  two  sets  of  material 
properties  were  assumed  to  the  represent  the  typical  behavior  of  composite  materials. 
Simulations  were  performed  on  laminates  with  both  sets  of  properties.  The  predicted  mechanical 
behavior  due  to  the  effect  of  oxidation  was  described.  This  included  illustrating  the  initiation  and 
progression  of  damage  in  the  laminate.  A  parametric  study  was  also  performed  to  determine  the 
effect  of  the  number  of  plies  on  the  mechanical  behavior  under  oxidation. 
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10.  Conclusions  and  Future  Work 

This  research  work  has  contributed  in  various  ways  to  help  develop  a  better  understanding  of 
textile  composites  and  materials  with  complex  microstructures  in  general.  An  instrumental  part 
of  this  work  was  the  development  of  a  framework  that  made  it  convenient  to  perform 
multiscale/multiphysics  analyses  of  advanced  materials  such  as  textile  composites  with  complex 
microstructures.  In  addition  to  the  studies  conducted  in  this  work,  this  framework  lays  the 
groundwork  for  continued  research  of  these  materials.  In  addition  to  implementing  an  oxidation 
model,  the  framework  was  also  used  to  implement  strategies  that  expedited  the  simulation  of 
oxidation  in  textile  composites  so  that  it  would  take  only  a  few  hours.  Finally,  a  coupled 
diffusion/oxidation  and  damage  progression  analysis  was  implemented  that  was  used  to  study  the 
mechanical  behavior  of  textile  composites  under  mechanical  loading  as  well  as  oxidation.  The 
following  sections  discuss  the  conclusions  drawn  from  the  work  performed  to  achieve  the 
objectives  of  this  research  effort.  This  section  concludes  by  mentioning  some  suggestions  for 
possible  future  work. 

10.1.  Simulation  of  oxidation  in  textile  composites 

The  oxidation  behavior  of  textile  composites  was  simulated  using  the  finite  element  framework 
that  was  developed  as  part  of  this  work.  This  involved  implementing  various  strategies  because 
of  the  multiple  scales  of  microstructure  involved  in  the  configuration.  An  oxidation  model  was 
implemented  based  on  the  model  developed  by  Pochiraju  et  al  to  simulate  oxidation  in  neat 
PMR-15  resin.  Homogenized  oxidation  material  properties  for  a  unidirectional  laminate  or  tow 
were  determined  assuming  that  the  fiber  was  impermeable  and  un-oxidizable.  The  homogenized 
properties  were  validated  using  different  configurations.  It  was  also  determined  that  the 
oxidation  behavior  in  heterogeneous  configurations  is  complex  and  depends  on  various  factors 
such  as  the  location  of  the  material  boundaries.  The  oxidation  model  had  severe  limitation  on  the 
element  size  and  time  step  size  based  on  the  finite  element  formulation.  Therefore,  a  typical 
oxidation  analysis  was  very  computationally  intensive  and  it  was  not  feasible  to  simulate 
oxidation  of  a  textile  composite  without  strategies  to  expedite  the  analysis.  Optimal  element  sizes 
were  determined  and  the  time  step  size  was  ramped  up  to  achieve  better  efficiencies.  An  adaptive 
meshing  strategy  was  also  developed  that  exploited  certain  characteristic  of  the  oxidation 
behavior  to  reduce  the  size  of  the  problem.  The  adaptive  meshing  strategy  was  able  to  give 
computational  time  savings  of  over  60%.  However,  these  initial  strategies  were  not  enough  to 
make  a  full  3D  oxidation  analysis  feasible.  Therefore,  a  decoupled  subdomain  strategy  was 
developed  that  divided  up  a  3D  analysis  domain  into  an  array  of  ID  domains  which  could  then 
be  solved  in  a  more  timely  manner.  The  ID  models  could  also  be  analyzed  independently  on 
different  processors  in  a  multi-core  machine  thereby  increasing  the  efficiency  even  further.  The 
decoupled  subdomain  strategy  was  validated  and  used  in  conjunction  with  the  adaptive  meshing 
strategy  to  simulate  oxidation  of  a  plain  weave  laminate.  The  analysis  revealed  that  the  tow 
architecture  of  the  textile  composite  had  a  significant  effect  on  the  distribution  of  the  oxidation 
front  location  in  the  woven  configuration.  After  200  hours  of  oxidation  of  a  200  micron  thick 
ply,  the  smallest  oxidation  layer  thickness  was  found  to  84  microns  whereas  the  largest  was  110 
microns.  Later  increases  in  efficiency  and  the  use  of  parallel  algorithms  eventually  led  to  the 
ability  to  perform  three-dimensional  oxidation  analysis  in  a  timely  manner.  However,  when  3D 
results  are  compared  to  results  obtained  through  the  decoupled  subdomain  strategy,  it  is  seen  that 
the  strategy  performs  remarkably  well  given  the  significant  simplifying  assumptions  involved. 
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10.2.  Prediction  of  damage  progression  in  textile  composites  under 
oxidation 

The  multiscale/multiphysics  analysis  framework  was  used  to  implement  a  coupled 
diffusion/oxidation  and  continuum  damage  analysis  to  study  the  mechanical  behavior  of  textile 
composites  in  oxidizing  environments.  The  current  implementation  of  the  coupled  model 
considers  only  the  effect  of  oxidation  on  the  mechanical  behavior  and  not  vice  versa.  Since  only 
one-way  coupling  was  assumed,  the  oxidation  simulation  could  be  performed  independently 
before  the  coupled  analysis.  The  coupled  analysis  was  used  to  predict  progressive  damage  in  a 
Graphite/PMR-15  plain  weave  laminate  that  is  uniaxially  loaded  to  a  fixed  strain  level  and  then 
exposed  to  oxidation  through  the  top  and  bottom  surfaces.  A  constitutive  model  was  developed 
that  degrades  the  engineering  properties  depending  on  the  mechanical  state  and  how  much  the 
material  has  oxidized.  Due  to  lack  of  a  full  set  of  reliable  material  properties,  two  sets  of  material 
properties  were  assumed  to  represent  the  typical  behavior  of  composite  materials.  The  predicted 
mechanical  behavior  due  to  the  effect  of  oxidation  was  described  and  an  attempt  was  made  to 
explain  some  of  the  behavior  observed.  The  simulations  showed  the  evolution  of  damage  in  the 
composite  as  it  undergoes  oxidation.  It  was  seen  that  in  some  cases  the  effect  of  oxidation  is  not 
localized  and  that  damage  also  occurs  in  regions  that  are  not  oxidized  due  to  load  redistribution. 
The  simulations  also  showed  the  variation  of  the  volume  averaged  stress  in  the  laminate  as  the 
laminate  oxidizes.  It  was  seen  that  the  various  stress  allowables  of  the  materials  in  the  laminate 
had  an  effect  on  this  behavior.  A  parametric  study  was  also  performed  to  determine  the  effect  of 
the  number  of  plies  on  the  mechanical  behavior  under  oxidation.  The  simulations  predicted  a 
significant  drop  in  the  load  carried  by  a  2-ply  laminate  for  different  strain  levels  and  the  load 
drop  reduced,  as  expected,  when  the  number  of  plies  was  increased  to  4  and  6.  However,  the 
proportion  by  which  the  load  drop  reduces  was  not  very  intuitive  and  indicates  that  the  material 
properties  and  the  microstructure  of  the  textile  laminates  have  a  complicated  effect  on  the 
behavior  under  oxidation. 

10.3.  Future  Work 

Over  the  course  of  this  research  work,  several  ideas  came  up  that  might  have  successfully  helped 
in  advancing  the  understanding  of  these  advanced  materials.  However,  not  all  of  them  could  be 
pursued  due  to  various  reasons.  In  addition  to  this,  there  are  some  obvious  extensions  to  the 
research  work  presented  in  this  dissertation.  Some  of  them  are  listed  below: 

1.  For  analyzing  even  smaller  length  scales,  hybrid  models  directly  linking  atomistic  regions 
to  continuum  finite  element  regions  have  been  developed  by  several  researchers.  These 
include  the  FEAt  model  [59],  the  MAAD  approach  [60,61],  the  QuasiContinuum  method 
[62-64]  and  the  coupled  atomistic  and  discrete  dislocation  plasticity  (CADD)  approach 
[65]  and  the  bridging  scale  method  [66].  Currently,  the  hierarchical  strategies  explained 
in  this  work  are  implemented  only  for  the  continuum  mechanics  regime.  However,  it 
might  be  worthwhile  to  explore  the  possibility  of  using  these  strategies  in  analyzing 
multiple  scale  problems  involving  the  atomistic  scale. 

2.  The  current  finite  element  formulation  for  the  oxidation  model  assumes  that  the  time  step 
is  small  enough  that  the  assumptions  to  account  for  the  nonlinearity  hold.  Future 
extensions  to  the  model  could  look  at  defining  a  residual  and  iterating  to  drive  the 
residual  to  zero  to  account  for  the  nonlinearity  at  each  time  step. 
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3.  When  trying  to  replace  a  heterogeneous  material  with  a  homogenized  material  in  an 
oxidation  model,  it  is  reasonable  to  assume  that  some  or  possibly  all  of  these  properties 
might  change.  In  addition  to  the  current  homogenization  strategy,  there  is  at  least  one 
other  possible  approach  for  achieving  this  goal.  One  is  to  use  a  multi-scale  analysis  that 
keeps  track  of  the  ‘local’  information  such  as  oxidation  state  and  actual  average 
concentration  in  the  constituent  matrix  in  the  homogenized  material.  Given  this 
information,  it  would  be  possible  to  calculate  the  reaction  rate  R  at  a  particular  material 
point  in  the  tow’s  constituent  matrix  using  Eq.(3.55).  Next,  the  effective  reaction  rate  for 
the  larger  scale  homogenized  tow  is  determined  by  a  simple  rule  of  mixtures  and  plugged 
into  the  governing  equations.  When  the  equations  for  a  time  step  are  solved,  the 
calculated  concentrations  are  transformed  back  to  the  local  scale  using  a  rule  of  mixtures 
in  order  to  keep  track  of  the  oxidation  state  of  the  constituent  matrix.  Thus,  a  continuous 
transfer  of  information  between  the  two  scales  needs  to  be  maintained  throughout  the 
simulation. 

4.  In  this  work,  the  fibers  in  the  tows  are  idealized  to  be  in  a  square  array  and  the  fibers  are 
assumed  to  be  impermeable  and  do  not  oxidize.  While  there  are  other  factors  that  can 
influence  the  oxidation  behavior  in  composites  such  as  the  properties  of  the  fiber/matrix 
interface  or  interphase,  they  are  not  taken  into  account  for  the  homogenization  model 
described  in  this  work.  Cracks  in  the  matrix  or  along  the  fiber/matrix  interface  can  also 
affect  the  oxidation  behavior  by  allowing  rapid  ingress.  Depending  on  the  type  of  damage 
that  is  observed  in  these  composites,  it  might  be  possible  to  account  for  their  effects  in 
the  homogenized  model.  For  example,  if  the  damage  is  diffuse,  the  homogenized 
properties  can  be  degraded  appropriately  or  if  the  damage  is  confined  to  certain  areas, 
cracks  can  be  introduced  in  the  homogenized  model.  These  and  other  such  factors  should 
be  addressed  in  a  likely  extension  to  the  model. 

5.  The  oxidation  level  information  from  the  decoupled  subdomain  strategy  is  currently 
approximated  as  a  single  linear  function  to  define  the  active  zone.  A  better  approximation 
could  be  made  using  a  few  more  points  to  define  a  piecewise  linear  function  for  the  active 
zone. 

6.  A  simple  constitutive  model  or  property  degradation  scheme  was  developed  to  account 
for  the  effect  of  oxidation  on  the  mechanical  behavior.  This  scheme  can  be  modified  and 
enhanced  in  the  future  when  the  effect  of  oxidation  on  the  coupled  oxidation-mechanical 
behavior  is  more  accurately  determined.  This  can  also  include  a  constitutive  model  to 
account  for  the  effect  of  mechanical  damage  on  the  oxidation  behavior,  which  would 
make  the  analysis  fully  coupled. 

7.  The  effects  of  thermal  expansion  and  the  shrinkage  of  the  matrix  due  to  oxidation  are 
ignored  in  the  current  implementation  of  the  coupled  model.  These  are  effects  that  should 
to  be  considered  in  future  implementations  of  the  coupled  analysis  model  in  order  to 
more  accurately  represent  the  behavior  of  the  underlying  mechanisms. 
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